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In  this  research,  we  develop  and  apply  correlated  methods  to  infinite,  periodic 
systems,  especially  polymers,  to  study  the  effects  of  electron  correlation  on  their 
properties. 

Second-order  MBPT  (MBPT(2))  is  implemented  for  infinite,  periodic  systems. 
The  convergence  of  the  MBPT(2)  contributions  to  the  total  energy  per  unit  cell 
and  band  energy  with  lattice  summations  and  other  numerical  cutoffs  are  assessed 
numerically. 

Using  our  developed  MBPT(2)  program,  the  correlated  quasi-particle  energies 
for  valence  bands  of  polyethylene  are  calculated.  The  MBPT(2)  band  energies  accu- 
rately explain  the  measured  photoelectron  spectra  (XPS,  UPS)  of  polyethylene  and 
resolve  long-standing  disagreements  among  three  experiments. 
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We  have  calculated  the  equilibrium  structures  and  vibrational  frequencies  of 
polymethineimine.  The  latter  represents  the  first  time  an  ab  initio  correlated  method 
has  been  applied  to  this  problem  in  an  infinite  system.  Our  numerical  vibrational 
frequencies  agree  with  measured  values  very  well. 

We  have  also  calculated  the  band  gaps,  the  most  important  property  of  semi- 
conducting materials,  for  polyacetylene  and  polymethineimine  using  MBPT(2).  We 
found  that  MBPT(2)  improves  the  Hartree-Fock  band  gaps  greatly  for  both  systems 
compared  to  experimental  observations. 

We  have  extended  many-body  perturbation  theory  (MBPT)  to  give  any  order 
corrections  to  IP's  and  EA's  for  both  finite  and  infinite  systems,  which  permits  a 
direct  evaluation  of  correlation  corrections.  We  show  that  only  two  types  of  dia- 
grams remain  after  the  mutual  cancellations  among  the  contributions  of  the  orbital 
relaxation  and  electron  correlation. 

We  show  rigorously  that  MBPT,  coupled  cluster  theory,  and  many-body  Greens 
function  or  propagator  methods  all  converge  uniformly  with  lattice  summations  in 
periodic  systems  including  polymers,  surfaces,  and  crystals.  Besides  the  importance 
in  establishing  the  correct  scaling  properties  for  these  methods  in  inhomogeneous 
systems,  the  proof  provides  useful  guidance  for  lattice  summation  cutoffs. 
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CHAPTER  1 
INTRODUCTION 

1.1    Importance  of  Electron  Correlation 

The  numerically  and  conceptually  correct  description  of  the  instantaneous  Cou- 
lombic  interactions  among  many  electrons,  the  correlation  problem,  has  been  the  focal 
point  of  atomic  and  molecular  theory  for  many  years.1  Enormous  progress  has  been 
made  in  the  treatment  of  electron  correlation  in  molecules,2  and  in  uniform  model 
extended  systems  like  the  electron  gas3  or  nuclear  matter.4  However,  the  correla- 
tion problem  in  real,  nonuniform  extended  systems  like  crystals,  polymers,  or  other 
solid  state  materials  remains  one  of  the  foremost  problems  in  the  field  of  electronic 
structure.5  In  these  real,  non-uniform  extended  systems,  it  is  the  electron  correlation 
that  is  critical  to  a  description  of  many  of  the  most  interesting  problems  such  as  the 
band  gaps,  optical  and  magnetic  properties,  superconductivity,  etc. 

Currently  the  most  widely  used  method  in  solid  state  applications  is  density 
function  theory  (DFT)  within  the  local  density  approximation  (LDA)6"8.  Since  the 
density  is  a  much  simpler  object  than  a  wavefunction  (or  the  density  matrix),  it  offers 
an  appealing  route  toward  implementation.  Generally  speaking,  as  implemented,  the 
DFT  method  approximately  includes  electron  correlation  and  gives  better  results  than 
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Hartree-Fock  (HF)  theory.  But  it  still  fails  for  many  problems,  notably  for  "band- 
gaps."9  Furthermore,  there  is  no  technique  to  systematically  improve  the  functionals 
used  in  DFT  and,  thus,  no  clear  path  to  converged  results. 

An  alternative  route  to  converged,  correlated  results  is  offered  by  Many-body 
perturbation  theory10  (MBPT)  and  its  infinite-order  generalization,  coupled  clus- 
ter (CC)  theory.11  Built  on  top  of  a  suitably  chosen  independent  particle  model, 
CC/MBPT  provides  a  systematic  way  to  obtain  the  essential  effects  of  correlation. 
This  most  accurate  results  for  molecules  are  frequently  of  this  type.12  Furthermore, 
because  of  their  correct  scaling  with  size  (i.e.  size  extensivity),  unlike  traditional  con- 
figuration interaction  (CI)  methods,  CC/MBPT  offer  appropriate  tools  for  extended 
systems.  Propagator13-18  or  Greens  function  methods  (GFM), 19-23  provide  another 
correlated  tool  to  calculate  the  electron  correlation  corrections  to  band  energies.  Now 
with  the  convergence  of  reliable  HF  packages24-32  for  extended  systems,  the  next  step 
toward  predictive  results  is  to  include  electron  correlation. 

1.2    Previous  Developments  in  Extended  Systems 

Early  efforts  to  account  for  electron  correlation  in  extended  systems  were  lim- 
ited to  the  development  of  various  phenomenological  models.  The  first  quantum- 
mechanical  based  model  was  provided  by  Toyozawa33  who  considered  the  interaction 
of  an  extra  particle  with  the  rest  of  the  electrons  by  introducing  the  concept  of  the 
electronic  polaron.  Toyozawa's  model  was  generalized  for  the  interaction  of  an  elec- 
tron and  a  hole  by  Haken  and  Schottky34  and  applied  to  alkali-halide  crystals  within 
the  framework  of  the  tight-binding  approximation  by  Inoue  et  al.35  Hermanson36 
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introduced  another  model  for  insulators  in  which  a  plasmon  field  represents  the  col- 
lective excitations  of  a  valence-electron  system. 

The  alternative  ab  initio  viewpoint  has  seldom  been  pursued,  even  though  the 
simplest  such  method  (MBPT(2))  has  been  shown  for  molecules  to  typically  account 
for  90%  of  the  correlation  energy,  and  to  offer  a  factor  of  two  improvement  in  other 
properties.12  Such  second-order  many-body  perturbation  theory  (MBPT(2))  correc- 
tions to  HF  have  been  considered  by  Kunz  and  coworkers,37  Suhai,38  and  Liegener.39 
Suhai40-42  has  done  MBPT(2)  calculations  for  structures  and  band  gaps  in  polyacety- 
lene  and  some  other  systems.  From  the  inverse  Dyson  equation  with  the  irreducible 
self-energy  part  in  the  diagonal  approximation,  Liegener39  has  calculated  the  second- 
and  third-order  corrections  to  the  band  gap  for  alternating  trans-polyacetylene.  How- 
ever, there  is  a  large  inconsistency  between  Suhai's40'41  and  Liegener's39  results  that 
requires  resolution. 

1.3    Implementation  of  MBPT(2)  and  Its  Applications  to  Polymers 

In  this  research,  we  would  like  to  develop  and  extend  the  correlated  methods  suc- 
cessfully used  in  finite  molecules  to  periodic,  infinite  systems  and  use  them  to  explain 
and  predict  the  properties  such  as  equilibrium  geometries,  band  gaps,  photoelectron 
spectra,  vibrational  frequencies,  of  polymers.  As  a  first  step,  we  have  implemented 
MBPT(2)  for  one  dimensional  periodic,  infinite  systems  and  applied  the  program  to 
polyacetylene,  polyethylene,  and  polymethineimine. 

In  chapter  2,  a  brief  description  of  the  HF  method  for  extended  systems  will 
be  given  and  then  the  detailed  formulas  for  MBPT(2)  corrections  to  both  the  total 
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energy  per  unit  cell  and  band  energies  will  be  derived.43  We  will  carefully  study  the 
convergence  behaviors  of  MBPT(2)  with  lattice  summations  cutoff  N,  the  number  of 
points  K  taken  in  the  integrals  over  reciprocal  space  and  cutoff  threshold  10_c  for 
two-electron  atomic  orbital  integrals.  We  will  show  that  the  MBPT(2)  correlation 
correction  to  the  band  structure  converges  very  slowly  with  N  and  demands  a  large  K 
while  the  MBPT(2)  correction  to  the  total  energy  per  unit  cell  converges  much  faster 
with  N  and  needs  a  much  smaller  K  while  neither  MBPT(2)  correction  is  sensitive  to 
the  cutoff  of  the  two-electron  atomic  orbital  integrals  when  C  >  5.  We  will  also  show 
that  the  previous  MBPT(2)  results  were  obtained  either  with  too  few  unit  cells  such 
that  the  convergence  with  N  had  not  been  reached,  or  that  the  zeroth-order  Hartree- 
Fock  results  were  inadequately  converged.  MBPT(2)  with  a  DZP  basis  improves 
the  Hartree-Fock  band  gap  from  5.57  eV  to  3.22  eV  at  the  experimentally  estimated 
geometry,  compared  to  the  measured  ~2  eV  peak  in  the  absorption  spectrum44  of 
the  system. 

X-ray  (XPS)45-48  and  ultraviolet  photoelectron  spectroscopy  (UPS)49-52  provide 
rich  information  about  the  valence  bands  of  extended  systems.  The  further  develop- 
ment of  angle-resolved  UPS  (ARUPS)  can  even  be  used  to  directly  observe  the  band 
structures.  In  the  last  two  decades,  these  two  methods  have  been  frequently  used 
to  elucidate  the  valence  electron  structures  of  synthetic  organic  polymers.  Among 
them,  polyethylene  is  well  studied.46-48'48-52  In  Chapter  3,  we  will  apply  MBPT(2)  to 
polyethylene.  We  will  show  that  electron  correlation  effects  for  the  valence  bands  in 
polyethylene  vary  from  1.5  eV  to  5.4  eV.53  The  correlated  quasi-particle  band  energies 
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given  by  MBPT(2)  with  a  polarized  6-31G**  basis  accurately  explain  the  measured 
photoelectron  spectra  (XPS,  UPS)  of  polyethylene  and  resolve  long-standing  disagree- 
ments among  these  experiments.  DFT  fails  to  provide  agreement  with  experiment. 
This  example  demonstrates  the  critical  role  of  correlated  ab  initio  theory  in  obtaining 
accurate  band  structures  for  extended  systems. 

Polymethineimine  is  a  conjugated  polymer,  isovalent  with  polyacetylene,  which 
was  first  synthesized  in  1971. 54  It  has  been  found  to  be  a  semiconductor.  Except 
for  its  vibrational  frequencies  which  have  been  measured,55  little  about  the  polymer 
is  known.  In  Chapter  4,  MBPT(2)  is  applied  to  polymethineimine  to  calculate  its 
equilibrium  structure,  band  gap,  and  vibrational  frequencies.56  The  latter  represents 
the  first  time  an  ab  initio  correlated  method  has  been  applied  to  this  problem  for 
an  infinite  system.  We  will  show  that  both  basis  set  and  electron  correlation  have  a 
strong  influence  on  its  optimized  geometry.  The  MBPT(2)  band  gap  with  a  polarized 
basis  set  (6-31G**)  is  4.7826  eV,  compared  to  8.54  for  Hartree-Fock.  Assuming  the 
same  error  as  for  polyacetylene,43  we  predict  a  gap  of  ~2.6  eV.  We  will  also  demon- 
strate that  electron  correlation  has  an  important  effect  on  the  computed  vibrational 
frequencies.  The  theoretical  vibrational  frequencies  calculated  with  the  6-31G  basis 
match  the  measured  data  very  well,  while  the  agreement  is  somewhat  poorer  in  a 
better  (6-31G**)  basis. 

1.4    Progress  in  Formal  Theory 

Besides  the  development  of  MBPT(2)  for  one-dimensional,  periodic,  infinite 
systems  and  the  application  of  the  method  to  polymers,  we  have  made  progress 
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in  formal  theory.  We  have  derived  a  general,  explicit  expression  for  quasi-particle 
energies  which  are  the  ionization  potentials  (IP)  for  occupied  orbitals  and  electron 
affinities  (EA)  for  unoccupied  orbitals,  respectively,  and  proved  that  MBPT,  CC, 
and  many-body  green  function  (MBGF)  methods  all  converge  uniformly  with  lattice 
summations,  although  the  integrand  for  the  integration  over  the  reciprocal  lattice 
vector,  k,  could  become  infinite  at  special  k  values. 

MBPT  has  been  well  developed  for  the  total  energy.  It  has  also  been  used  to 
calculate  ionization  energies  and  electron  affinities  but  limited  to  second-  and  third- 
order.  In  chapter  5,  we  will  extend  the  theory  to  give  any  order  corrections  to  IP's  and 
EA's  for  finite  and  infinite  systems,  which  permits  a  direct  evaluation  of  correlation 
corrections.  We  will  show  that  the  diagrams  for  IP's  and  EA's  can  be  classified 
into  two  types.  The  first  type  of  diagrams  are  linked  diagrams  which  can  be  further 
separated  into  two  sets.  One  set  consists  of  the  diagrams  having  one  and  only  one 
bubble  with  a  fixed  index.  They  mainly  account  for  orbital  relaxation  effects  that  are 
shown  to  have  a  finite  contribution  in  infinite  periodic  systems.  The  other  diagrams 
are  derived  from  the  diagrams  for  the  total  energy.  In  each  of  those  diagrams,  there 
is  one  and  only  one  line  with  a  fixed  index.  The  second  type  are  the  diagrams  which 
have  two  or  more  lines  with  fixed  indices.  Because  of  mutual  cancellations  among  the 
contributions  of  the  orbital  relaxation  and  electron  correlation,  the  only  remaining 
diagrams  are  those  which  fall  apart  by  switching  any  two  ending  vertices  of  the 
lines  that  have  fixed  indices.  The  connection  between  the  first  type  of  diagrams  of 
the  MBPT  correction  and  the  total  self  energy  with  the  diagonal  approximation  for 
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E  —  ej,0'  in  propagator13-18  or  Greens  function  methods19-23  is  explicitly  established 
for  any  order.  We  will  also  show  that  the  orbital  relaxation  has  a  finite  contribution  to 
the  band  energy  for  infinite  periodic  systems,  although  its  effect  on  the  zeroth-order 
orbital  energy  and  total  energy  per  unit  cell  is  infinitesimal. 

Recently,  the  possible  divergence  of  many-body  perturbation  theory  or  many- 
body  green  function  (propagator)  methods  with  lattice  summations  in  extended  sys- 
tems has  been  raised.58  The  convergence  of  these  methods  with  lattice  summations 
is  not  only  the  key  to  establishing  their  correct  scaling  properties59  in  inhomoge- 
neous  systems,  but  is  also  a  necessity  if  numerical  calculations  are  to  be  mean- 
ingful. In  Chapter  6,  we  will  rigorously  show  that  MBPT,10  CC,11-12  and  MBGF 
methods13-15-19-20  all  converge  uniformly  with  lattice  summations.60  Our  proof  is  given 
not  only  for  infinite  polymers  but  also  for  crystals.  In  practical  numerical  calcula- 
tions, one  should  not  pursue  the  convergence  of  the  integrand,  but  the  convergence 
of  the  correction  considered  with  lattice  summations. 


CHAPTER  2 

SECOND-ORDER  MBPT  FOR  EXTENDED  SYSTEMS  AND  ITS 
APPLICATION  TO  POLYACETYLENE 


2.1  Introduction 

In  this  chapter,  we  present  all  the  formulas  needed  to  implement  MBPT(2)  for 
infinite,  periodic  systems.  We  have  used  them  in  writing  our  MBPT(2)  programs 
for  polymers.  We  also  provide  an  in-depth  study  of  applying  MBPT(2)  to  polymers, 
with  application  to  trans-polyacetylene.  For  extended  systems,  the  calculated  results 
can  be  seriously  influenced  by  several  parameters  such  as  the  number  of  unit  cells 
(N)  taken  in  the  lattice  summations,  the  number  of  k-points  (K)  taken  in  the  first 
Brillouin  zone  in  the  integrals  over  k  and  the  cutoff  threshold  (10~c)  for  the  two- 
electron  atomic  orbital  integrals.  Suhai40  has  studied  the  convergence  of  the  second- 
order  correction  to  the  total  energy  per  unit  cell  (E$)  with  N.  However,  there  is  no 
serious  study  on  the  convergence  of  the  second-order  MBPT  correction  to  the  band 
structure  with  N.  There  is  also  no  study  about  the  proper  values  of  K  in  the  MBPT(2) 
total  energy  and  band  structure  calculations,  respectively,  and  the  influence  of  C  on 
the  numerical  results.  In  this  chapter,  we  will  show  that  the  MBPT(2)  correction 
to  the  total  energy  of  one  unit  cell  E$  converges  very  fast  with  N  and  is  not  very 
sensitive  to  either  K  or  C.  In  fact,  E<$  converges  much  faster  if  the  converged  HF 
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band  energies  and  wavefunctions  are  used  and  kept  unchanged,  than  in  Suhai's40 
calculation  where  he  took  N  in  the  HF  calculation  and  in  the  calculation  of  E$ 


structure  converges  very  slowly  with  N  and  demands  much  larger  K  values  to  give 
reliable  results.  We  provide  MBPT(2)  results  for  the  total  energy  and  band  gap  for 
alternating  trans-polyacetylene  at  several  geometries  and  basis  sets. 

The  plan  of  the  chapter  is  as  follows.  In  section  2.2,  the  HF  formulas  are 
summarized  to  facilitate  some  definitions.  In  section  2.3  and  2.4,  explicit  expressions 
for  the  MBPT(2)  total  energy  per  unit  cell  and  for  the  MBPT(2)  correction  to  the 
band  gap  are  given.  These  expressions  can  be  directly  used  for  implementation  and 
are  given  in  more  detail  than  previously.5'37'38  Implementing  these  expressions,  in 
section  2.5,  the  convergence  of  both  the  MBPT(2)  correction  to  the  total  energy  per 
unit  cell  and  that  to  the  band  energies  with  N,  K  and  C  are  studied.  In  section  2.6, 
we  apply  MBPT(2)  to  alternating  trans-polyacetylene.  We  close  with  concluding 
remarks  in  Sec.  2.7. 


Let  ai ,  a2  and  a3  be  the  three  basis  vectors  of  a  crystal  or  a  three-dimensional 
periodic  polymer.  Then  the  nonrelativistic  Hamiltonian  of  the  system  in  the  Born- 
Oppenheimer  approximation  is 


to  be  the  same.  Contrary  to         we  find  the  second-order  correction  to  the  band 


2.2    Hartree-Fock  Method 


H  = 


IRj  +  Rfl-Ri-R^I 


i  R.,,-4  Kj,B 

ZA 


(2.1) 
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where  =  iiai  +  i2a2  +  ^3  and  Rj  =  ji&i  +  j2a2  +  j3&s  are  the  lattice  vectors, 
Za  and  Zb  are  the  charges  of  the  nuclei  A  and  B  at  position  Ra  and  Rb  relative  to 
the  origins  of  their  unit  cells,  respectively,  and  the  prime  means  that  the  cases  where 
R,  =  Rj  and  A  =  B  simultaneously,  are  excluded. 

With  the  periodicity  of  the  Hamiltonian,  one  can  introduce  a  good  quantum 
number  k  =  k\bi  +  k2b2  +  k$bz  where  bj  (j  =  1, 2, 3)  are  determined  by  •  bj  =  5ij, 
i  =  1, 2, 3.  When  ki  (i  =  1, 2, 3)  are  limited  to  the  region  [— tt,  7r],  the  k  space  is  called 
the  first  Brillouin  zone  and  its  volume  is 

W   =   (27r)3b!  x  b2  •  b3  =  (27r)3/a!  x  a2  •  a3.  (2.2) 

According  to  the  Block  theorem,  the  one-electron  spatial  wavefunctions  can  be 
written  as 

^nk(r)   =   etk'runk(r),  (2.3) 

where  unk(r)  is  translationally  invariant.  ^nk(r)  can  also  be  expressed  as  the  sum- 
mation of  the  atomic  orbitals 


(2.4) 
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where 

A^r-R,)   =   EcpVP(r)  (2.5) 
p 

and 

4(t)   =   Xp(r-R,-Rp)  (2.6) 

are  the  atomic  orbitals  centered  at  Rp  in  the  jth  unit  cell.  N  in  Eq.(2.4)  is  the  total 
number  of  the  unit  cells  in  the  system  and  it  approaches  +oo.  The  overlap  between 
the  molecular  orbitals  is  given  by 

(Mr)\Mr))    =  E(^k)^r^EE^'Rj-,kR,(xP(r)|xJ9(r)) 

pq  R<  R.j 

=      £  (cfrcyX,  (2-7) 

pq 

where 

5pk9   =   E^kRj(xS(r)|x9(r))  (2.8) 
is  the  overlap  matrix.  From  the  orthonormalization  of  the  orbitals,  one  can  obtain 

E(C;krcf<Spk9  =  6nn,  (2.9) 

pq 
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The  Hartree-Fock  equation  for  the  system  is 


5,32 


FkCnk     =    eHFskGnk)  (2.10) 


where 


3 

is  the  Fock  matrix.  In  Eq.(2.11),  T£  ,  V3q,  J3pq  and  K3vq  are  the  kinetic,  electron- 
nucleus,  Coulomb  and  exchange  energy  matrices,  respectively.  They  are  given  by32 

TU   =   Jx0P(r){-\v2]xl(r)dr,  (2.12) 


pi 

Kh,A 


Jpq 


hlrs 


(2.14) 


K3 

PQ 


hlrs 


{pr 


(2.15) 
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where 

W\rhsl)  =  j  j  XQpWq{v)\v  -  rrW)x's(rVrrfr'  (2.16) 

are  the  two-electron  integrals  and 

1  occ'  r  * 
Dm   =  dk(°pk)  C9nkeik(R'-R")  (2.17) 

is  the  electron  density  matrix. 

From  the  HF  Block  orbitals  0„k(r),  the  HF  ground  state  of  the  system  can  be 
constructed  as 


=   |...(nkw)...),  (2.18) 

where  u  is  the  spin  quantum  number  of  the  spin  orbitals.  The  HF  total  energy  per 
unit  cell  is  then 


=   E  [2^  +  2^  +  2^-/^]^ 


^  1  Zj^Z[j 


where  the  prime  means  that  the  cases  where  |Rj  -       +  RB|  =  0  are  excluded. 

In  Eqs.(2.8),  (2.11)  and  (2.19),  there  are  infinite  summations  in  the  overlap  ma- 
trix Sk,  the  Fock  matrix  Fk  and  the  total  energy  per  unit  cell  E"F .  There  are  also 
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infinite  summations  in  7^  ,  Vjq,  J^q,  and  K3pq.  In  real  calculations,  it  is  not  feasible  to 
do  an  infinite  summation  and  the  lattice  summation  has  to  be  truncated  somewhere, 
so  one  has  to  impose  a  cutoff,  somehow.  Fortunately,  the  infinite  lattice  summations 
in  Sk,  TJ  and  converge.61  When  we  combine  V,  JJ  and  the  internuclear  inter- 
actions together,  the  infinite  lattice  summations  are  also  convergent.62  Thus  these 
summations  can  be  truncated.  However,  different  periodic  HF  calculations  may  use 
different  truncation  strategies  and  thus  give  different  results  even  when  using  the 
same  basis  set  and  the  same  geometrical  configuration.32  If  the  truncation  has  been 
done  properly,  the  difference  between  the  results  obtained  by  any  two  such  programs 
should  be  reasonably  small. 

The  orbital  energies  e^kF  form  the  well-known  HF  energy  bands,  where  n  is  the 
band  index.  According  to  Koopmans'  theorem,63  the  band  energies  for  the  occupied 
bands  n  and  unoccupied  bands  n',  are  equal  to  the  corresponding  ionization  potentials 
and  electron  affinities,  respectively,  i.e. 

dF   =  E$  -  E<HNi%k),  (2.20) 

41   =  E^\n'k')-E^l  (2-21) 

if  the  orbital  wavefunctions  do  not  change  appreciably  during  the  removal  or  addition 
of  an  electron  to  the  system.  Since  there  is  an  infinite  number  of  electrons  in  infinite 
systems,  removing  or  adding  one  cannot  change  anything,  so  the  above  Koopmans' 
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assumption  is  valid  for  the  Hartree-Fock  bands  for  extended  systems.  This  differs 
from  the  molecular  case. 

When  (n'k')  and  (nk)  are  the  lowest  unoccupied  and  the  highest  occupied  or- 
bitals,  respectively,  AcWF(nk  -»  n'k')  reaches  its  minimum.  This  minimum  value  is 


For  a  one-dimensional  periodic  system  or  a  polymer  with  periodicity  a  in  the  z 
direction  i.e.  Rt  =  iaez,  k  becomes  kez  or  simply  a  scalar  k  and  its  first  Brillouin  zone 
volume  is  2-n/a.  With  these  two  replacements,  all  the  formulas  given  in  this  section 
and  in  the  following  sections  for  MBPT(2)  can  be  directly  applied  to  one-dimensional 
periodic  systems. 


From  perturbation  theory,  we  know  that  the  exact  total  energy  of  the  system 
can  be  expressed  as 


called  the  HF  band  gap  and  is  denoted  as  Eg  . 


2.3    MBPT(2)  for  the  Total  Energy  per  Unit  Cell 


E   =   £(°>  +  E™  +  £(2)  +  .... 


(2.22) 


Since 


EHF     =  E(0)+E(l)j 


(2.23) 


Eq.(2.25)  can  be  rewritten  as 


E   =   EHF  +  E®  +  .... 


(2.24) 
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According  to  Raleigh-Schrodinger  perturbation  theory,5,10 


Fm         1  V    \(IJ\ru(l-P)\AB)\2  f2  25, 

^  —      A  eHF    i   fHF         fHF         fHF  ' 

4  I Jab  e/       eJ       e^  fcs 


where  7,  J  and  ^4,  B  denote  the  occupied  and  unoccupied  spin  orbitals,  respectively, 
in  the  HF  ground  state  of  the  system.  Let's  use  /,  J,  A  and  B  to  denote  the 
corresponding  spatial  orbitals  of  7,  J,  A  and  B,  respectively,  then  for  a  closed  shell 
system,  Eq.(2.25)  can  be  rewritten  as 


£(2)  y  2\(IJ\r^\AB)\2  -  Re[(IJ\r^\AB)(BA\r^\IJ)] 

I  JAB  eI         eJ     ~  tA     ~  eB 


For  an  extended  system,  I  is  a  composite  index  (ikj),  where  kj  is  confined  to 
the  first  Brillouin  zone,  and  likewise  for  J,  A  and  B.  Using  the  wavefunctions  given 
by  Eq.(2.4)  for  /,  J,  A  and  B,  we  obtain 


(IJ\r^\AB)   =   N~2         £  e[-i(k,R.+kJRJ-ka.R0-k6.R6)] 

x(Ak-(ri  -  R)Ak'(r2  -  R^r^A^  -  Ra) 
xAk"(r2-R6)).  (2.27) 


By  introducing  the  transformations 

rm  =  rm-R,  m  =  1,2,  (2.28) 
R;   =   Rp-R,   p  =  j,a,b,  (2.29) 
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Eq.(2.27)  can  be  simplified  as 

(IJ\r^\AB)   =   i-^eipf-iCki  +  kj-ka-kfcJ.Ri)] 

7V  R, 

>4  £  exP[ -  i(k,- ■  r;  - k.  •  b;  - k*  •  r;)] 

x(AIk'(r'1)A^(r^  -  R})|rS|^W  -  K) 
xA^-R'J).  (2.30) 

The  summation  jj  £R.  exP[  ~  i(ki  +  k,  —  ka  -  k(,)  •  Rj)]  gives  a  non-vanishing  result 
only  if 

(ki  +  kj-ka-kft)   =  2(Zibi  +  l2b2  +  /3b3)7r  (2.31) 

for  any  given  integers  l\,  Z2  and  /3  (k  =  0,  ±1,  ±2, %  —  1,  2,  3)  e.g. 

eo;p[  -  i(ki  +  kj  -  ka  -  k6)  •  R,-]   =   1  (2.32) 

for  every  R.  Let's  define  T(x)  as  a  function  which  moves  the  variable  x  back  to  the 
first  Brillouin  zone  by  n-fold  translations  of  27rbi,  27rb2  and  27rb3  if  it  is  out  of  the 
region.  Then  for  any  given  ki5  ka  and  kb,  there  is  always  one  and  only  one 


kj   =  T(ka  +  k6-ki) 


(2.33) 
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in  the  first  Brillouin  zone,  which  satisfies  Eq.(2.30).  Therefore,  Eq.(2.29)  can  be 
rewritten  as 

(IJ\r^\AB)   =   (Jk.,T(ka+ki_kj)g(ija6kikak6)/iVJ  (2.34) 

where 


Q(ijabkikakb)   =   E  E  exp{  -  z[(ka  +  kb  -  k*)  •  Rj  -  ka  •  Ra  -  k6  •  R6]} 

x(Alk'(r1)Aj(k"+k6-kl)(r2  -  R^r^fa  -  Ra)Ak*(r2  -  R6)) 
=     Y,    exp[i{ka  •  Rj  -  (ka  +  k6  -  ki)  ■  R2  +  k6  •  R3)] 

Ri,R2,R3 

x  E  (C^YiC^^-^yC^C^ipq111  |rR2sRs).  (2.35) 

pqrs 


Substituting  Eq.(2.34)  into  Eq.(2.26),  we  get 


Ei2)   =   ^EE  {2\Q(ijabktkakb)\2  -  RelQiijabkiKk,,) 

ijab  k,  k„  k(i 

xQ-iijbdkikfa)]}/  [e£f  +  e^ka+k6_ki)  -  eZ  ~  4f]  •  (2.36) 


Since  the  vector  k  has  in  total  N  points  in  the  first  Brillouin  zone,  each  k  point 
has  a  volume  W/N,  where  W  is  the  volume  of  the  first  Brillouin  zone.  Then  the 
summation  £k  in  Eq.(2.36)  can  be  replaced  by  the  integration  (N/W)Jdk.  After 
those  replacements,  Eq.(2.36)  becomes 


e(2)  =  w^L^L^Ibz^2^^^^ 


19 


-Re[Q{ijabkikakb)Q*(ijbakikbka)}} 

I  (Cikf  +  ^ka+kj-ki)  _  €aka  ~  ebYb)  >  (2-37) 

where  BZ  means  the  integral  range  is  the  first  Brillouin  zone.  In  an  extended  system, 
the  total  energy  of  the  system  is  always  infinite  and,  therefore,  only  the  total  energy 
per  unit  cell  is  meaningful.  The  total  energy  per  unit  cell  is  defined  as 


EUC    —     Jr,  (2.38) 


which  can  also  be  expressed  as 


Euc  =  E™  +  £g>  +  ....  (2.39) 


From  Eqs.(2.37)  and  (2.38),  we  know  that 


b)\2 


Eu}   =   WlT,  [   dki[    d*a[  dkJ^Qiijabk.Kk, 

W    ijabJBZ         JBZ  JBZ 

-Re[Q(ijabkikakb)  x  Q*(tj6akik6kfl)]} 

/ (Cikf  +  ej7Xka+k6-k,)  _  eaka  ~  e6kf )  ■  (2.40) 

2.4    MBPTY2)  for  the  Band  Structure 
Beyond  the  HF  approximation,  one  can  still  define  the  band  energies  as37 

enk   =   E^-E^N-l\nk)  (2.41) 
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for  occupied  orbitals  and  as 


en,w    =   E^N+1\n'k')  -  E^N)  (2.42) 

for  unoccupied  orbitals.  The  enk  and  en>^>  are  no  longer  one-electron  energies,  but  are 
quasi-particle  energies.  The  HF  one-electron  energy  bands  are  a  first  approximation 
to  them.  One  can  also  define  the  energy  band  gap  Eg  for  the  quasi-particle  energy 
bands  as  the  energy  difference  between  the  lowest  unoccupied  orbital  energy,  en/k< ,  and 
the  highest  occupied  orbital  energy,  c„k.  Similarly,  the  HF  energy  band  gap  E^F  is  a 
first  approximation  to  the  quasi-particle  energy  band  gap,  Eg. 

Just  as  for  the  total  energy,  we  can  use  perturbation  theory  to  express  ep  as 

enk   =   e£)+e£)  +  e£)  +  ...,  (2.43) 


where  P  denotes  either  nk  or  n'k'.  From  Koopmans'  theorem,63  we  know  that 


e 


(0)     _  rHF 


P 


e"P",  (2.44) 


s(p]    =   0.  (2.45) 


Koopmans'  theorem  demands  that  the  wavefunction  of  the  N+l  electron  system  and 
that  of  the  N-l  electron  system  be  single  determinants,  respectively,  and  assumes 
that  the  orbital  relaxation  due  to  the  removal  or  addition  of  an  electron  to  the  N 
electron  system  be  negligible.  This  means  that  at  zeroth-order,  there  are  no  other 
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states  with  the  same  symmetry  degenerate  to  $$pl\n'k'u}')  =  (n'k'u/)|0)  and  to 
^^(■nku))  —  (riku)\0).  Let's  use  e"F  and  efff  to  denote  the  lowest  unoccupied  and 
highest  occupied  HF  orbital  energies,  respectively.  When  c^f  <  e"F  —  E"F ,  it  is  most 
likely  that  there  is  another  state  whose  symmetry  is  the  same  as  that  of  (nkcj)|0)  and 
which  is  degenerate  with  (nkcj)|0).  Then  non-degenerate  perturbation  theory  breaks 
down.  Similarly,  we  can  conclude  that  Eq.(2.43)  becomes  a  poor  approximation  for 
en'k'  when  e%F,  >  e^F  +  EffF.  Therefore,  in  our  following  discussions,  we  confine 
ourself  to  the  energy  bands  which  satisfy  the  following  condition 


elu    +  Lg      >  tP     >  eho    -  hg     •  \ZAb) 

Using  MBPT(2)  for  E^'^ink)  and  E^N+l\n'k'),  respectively,  in  a  com- 
mon set  of  orbitals  for  the  E^  system,  we  can  explicitly  subtract  .E^-1)  and  E^N+1^> 
to  give  a  formula  for  a  direct  evaluation  of  the  band-gap  instead  of  obtaining  it  as  a 
difference  of  two  large  numbers.64  Then  using  Eqs.(2.26)  and  (2.41),  we  can  get  the 
following  expression  for  the  MBPT(2)  correction  to  the  band  energies  e^37,57 

eg}   =   U(P)  +  V{P),  (2.47) 

where  P  can  be  either  an  occupied  or  an  unoccupied  orbital, 


U(P)        y     1\(PI\\AB)\2  -  Re[(PI\\AB)(BA\PI)} 

u^  >  2^,2^,  HF  ,  fHF      Thf      Thf  '  (^.48) 
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v(p)  _        ^i^i^P-Myii/^/iP,)!  (24g) 


Substituting  Eq.(2.34)  into  Eqs.(2.48)  and  (2.49),  we  obtain 

U(P)   =    (^)2£/    dka  [  dk6{2|Q(pza&kpkak6)|2 

^W  '     iabJBZ  JBZ 

-Re[Q(piabkpkakb)Q*  (pibakpkbka)]  ] 

/  (*hf  ,  MF  mf     ^hf\  /0  r^^ 

/  ^cpkp  +  ejT(kQ+k(,-k,)  _  eaka  _  c6k6  J  i  l^.DUj 


K(P)   =  ^2Y,JBzdkiJBzdk3{2\Q(paijkpklkj)f 


i  V 

-Re[Q(paijkpkikj)Q*  (pajikpkjki)]  } 


/  f  HF  I     HF  HF      CHF\  /Q  r,\ 

/  ^pkP  +  ^TCki+k.-kp)  -  eik,   -  cjk,  J  •  {t-Ol) 


Both  C/(P)  and  V(P)  can  be  rewritten  in  more  symmetric  forms 


U{P)   =    (w)^f    dK  I  dkb{\Q(piabkpkakb)\2 

^  W  '     iab    BZ  JBZ 

+\Q(pibakpkbka)\2  -  Re[Q(piabkpkakb)Q*(pibakpkbka)]} 

I  (epkp  +  ejroca+kt-ki)  ~  e"ka  ~  e"kb)  i  (2.52) 

v{p)  =  (i)2^/B/k'/fl/k4i^^'k^)i2 

Q.IJ 

+  \Q{pajikpkJki)\2  -  ReiQipaijkpkikjWipajikpkjki)]} 

/  (cpkp  +  ^ki+kj-kp)  -       -  eg*")  ■  (2.53) 
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From  Eqs.(2.52)  and  (2.53),  we  see  that  all  numerators  in  both  equations  are  always 
positive.  Under  the  condition  given  in  Eq.(2.46)  for  epkp,  all  denominators  in  V(P) 
are  positive  while  those  in  U(P)  are  negative.  Therefore,  with  the  same  conditions, 
V(P)  is  always  positive  and  U (P)  is  always  negative,  that  is, 


When  P  is  an  occupied  orbital,  it  is  easy  to  see  that  the  denominators  of  V(P)  in  most 
cases  are  smaller  than  the  absolute  values  of  £/(P)'s  denominators.  If  numerators 
V(P)  and  U(P)  are  comparable,  V(P)  should  be  larger  than  the  absolute  value 
of  U(P).  Then,  the  second-order  correction  e„2^  for  valence  bands  are  positive  or 
the  correction  moves  bands  upward.  Similarly,  we  can  see  that  the  second-order 
correction  e^,  for  conduction  bands  are  negative,  the  correction  moves  the  bands 
downward.  Because  of  this,  the  second-order  MBPT  energy  band  gap  is  smaller  than 
the  HF  energy  band  gap. 

The  Green's  function  method  (GFM)  provides  another  approach  to  calculate 
the  quasi-particle  energy  bands.  The  Dyson  equation  provides  the  exact  E(N  ±  1) 
energies  in  a  formally  one-particle  picture.4  From  the  inverse  Dyson  equation,  with 
the  irreducible  self-energy  part  in  the  diagonal  approximation  being  correct  to  second- 
order,  one  can  get 


U(P)  <  0, 


V(P)  >  0. 


(2.54) 


EE 


I J  A 


2\(PA\\IJ)\*  -  Re[(PA\\IJ){JI\PA)] 
ef  +  e?F  +  e§F-e?F-e?F 
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+  EE 


2\{PI\\AB)\2  -  Re[(PI\\AB)(BA\PI)} 


(2.55) 


/  AB 


Comparing  Eq.(2.55)  with  Eqs.(2.47),  (2.48)  and  (2.49),  we  see  that  the  only  dif- 
ference between  the  two  formulas  is  that  Eq.(2.55)  has  an  extra  term  eft  in  its 
denominators.  As  an  iterative  solution  of  Eq.(2.55),  the  first  iteration  corresponds  to 
Eqs.  (2.47-2.49).  Higher  iterations  introduce  some  higher-order  correlation  effects. 


Following  the  formulas  given  in  section  2.4,  we  have  written  a  second-order 
MBPT  program,  PLMBPT,  for  extended  systems.  In  this  section,  we  will  study  how 
the  MBPT(2)  correction  to  the  total  energy  per  unit  cell,  E$,  and  the  MBPT(2) 
correction  to  the  band  structures,  e^,  vary  with  the  number  of  unit  cells,  N,  in 
the  lattice  summations,  cutoff  of  the  two-electron  integrals,  10~c,  and  number  of 
k-points,65  K,  taken  for  the  integral  over  k  in  the  first  Brillouin  zone.  The  objective 
is  to  discuss  the  dependence  of  the  calculated  E$  and  e^l  to  establish  reliable  values 
at  a  given  precision. 

There  are  several  HF  packages  for  extended  systems.24-32  In  Appendix  A,  three 
of  them  have  been  compared.  In  our  following  discussion,  we  choose  PLH93,  devel- 
oped by  the  Andre's  group,32  to  provide  the  HF  wavefunctions  and  the  HF  orbital 
energies  for  our  MBPT(2)  and  GFM  calculations.  We  will  fix  N  to  be  21,  e.g.  10 
unit  cells  at  each  side  of  the  reference  cell,  and  the  threshold  of  SCF  convergence  to 
be  10"8  in  our  HF  calculations  except  where  indicated. 


2.5    Numerical  Convergence  of  MBPT(2) 
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Table  2.1.  Structural  parameters  of  the  alternating  trans-polyacetylene  chain  (dis- 
tances in  A,  angles  in  degree). 


ri 

T2 

r3 

a 

7 

Gla 

1.3276 

1.4714 

1.0833 

123.97 

118.25 

G2a 

1.3719 

1.4824 

1.1026 

123.64 

118.34 

G3b 

1.3490 

1.4610 

1.0900 

123.80 

118.10 

G4b 

1.3650 

1.4500 

1.0900 

124.00 

118.05 

G5b 

1.3660 

1.4500 

1.0900 

123.90 

118.05 

G6a 

1.3634 

1.4450 

1.0872 

123.75 

118.13 

G7C 

1.3800 

1.4300 

1.0900 

122.00 

118.50 

a  Taken  from  reference  41. 
b  Taken  from  reference  40. 
c  Taken  from  reference  66. 


We  will  use  alternating  trans-polyacetylene  as  the  extended  system  in  our  dis- 
cussion. The  structure  and  the  structural  parameters  r1}  r2,  r3,  a  and  7  used  here  are 
the  same  as  those  used  by  Suhai.41  Table  2.1  summarizes  the  geometries  used  in  this 
section  and  the  following.  Geometries  from  Gl  to  G6  in  Table  2.1  were  optimized  by 
Suhai  with  different  basis  sets.40'41  G7  was  an  experimentally  estimated  geometry.66 
To  study  convergence,  we  use  G2  as  the  geometry  for  the  system  and  STO-3G  as  the 
basis  set  in  this  section.  Then  we  will  extend  the  basis  to  DZP  for  a  serious  study  in 
Section  VI. 

2.5.1    Numerical  convergence  of  MBPT(2)  for  the  total  energy  per  unit  cell 

In  this  subsection,  we  will  concentrate  on  the  convergence  of  E$  with  N,  C 
and  K. 
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Table  2.2.  Variation  of  E$  with  the  lattice  summation  cutoff  N.  System:  alternating 
trans-polyacetylene;  basis  set:  STO-3G;  geometry:  G2;  energy  unit:  atomic  unit.  For 
PLH93,  N=21.  For  PLMBPT,  C=6  and  K=21.  The  two  Is  core  bands  have  been 
omitted. 


N=5 

N=9 

N=15 

-c'uc 

-0.12344138 

-0.12273885 

-0.12252190 

Table  2.2  shows  how  the  correlation  energy,  E$,  of  alternating  trans  polyacety- 
lene  varies  with  N.  The  E@)  obtained  with  N  =  5  is  already  about  99.43%  of  that 
obtained  with  N  =  15.  When  N  —  9,  the  percentage  reaches  99.8%.  Tims,  we  can 
conclude  that  E^J  convergences  rapidly  with  N.  Suhai  made  a  similar  study.40  He 
took  poly-NH  as  the  extended  system.  In  his  calculations,  he  let  N  be  the  same  in 
both  the  HF  and  MBPT(2)  calculations,  i.e.  N  in  HF  also  changes.  Then  E$  in  his 
calculations  did  not  converge  as  fast  as  E$  does  in  ours.  In  fact,  he  obtained  less 
than  94%  of  the  total  E$  when  TV  =  10  (equivalent  to  N=5  in  our  case  since  our 
unit  cell  is  C2H2). 

Table  2.3  demonstrates  how  E"F  resists  variation  with  C.  In  fact,  all  three  E^ 
obtained  with  C—o.  6  and  7,  respectively,  are  consistent  in  the  first  six  significant 
figures.  Therefore,  it  is  adequate  to  take  C=5  in  the  calculation  of  the  total  energy 
per  unit  cell. 

Table  2.4  shows  how  E$  varies  with  K.  The  difference  between  E$  calculated 
with  K  =  11  and  K  =  21  is  about  0.2%  while  the  values  of  E$  obtained  with 
K  =  21  and  K  =  41  are  the  same  to  the  first  six  significant  figures.  This  means  that 
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Table  2.3.  Variation  of  E$  with  the  cutoff  10"c  of  two-electron  integrals.  Sys- 
tem: alternating  trans-polyacetylene;  basis  set:  STO-3G;  geometry:  G2;  energy  unit: 
atomic  unit.  For  PLH93,  N=21.  For  PLMBPT,  N=15  and  K=21.  The  two  Is  core 
bands  have  been  omitted. 


C=5 

C=6 

C=7 

P(2) 
•C'uc 

-0.12252173 

-0.12252190 

-0.12252160 

Table  2.4.  Variation  of  Eu2}  with  the  number  of  k-points  K  taken  in  the  integrals 
over  k  in  the  first  Brillouin  zone.  System:  alternating  trans-polyacetylene;  basis  set: 
STO-3G;  geometry:  G2;  energy  unit:  atomic  unit.  For  PLH93,  N=21.  For  PLMBPT, 
N=15  and  C=6.  The  two  Is  core  bands  have  been  omitted  in  the  summations  of  E$. 


K  =  ll 

K=21 

K=41 

E$  -0.12270513 

-0.12252190 

-0.12252166 

it  is  sufficient  to  take  K  =  21  in  very  accurate  E$  calculations.  In  fact,  K  =  11 
already  satisfies  the  precision  of  most  current  applications. 

2.5.2    Numerical  convergence  of  MBPT(2)  for  the  band  structure  and  band  gap 

In  this  subsection,  we  study  how  eft,  6p)  and  then  E^\  E^  vary  with  N,  C 
and  K.  We  will  show  that  they  have  much  slower  convergence  with  N  and  K  than 
does  the  total  energy  per  unit  cell.  In  MBPT(2)  calculations  of  the  band  structure,  N 
and  K  should  be  large  enough  to  guarantee  that  the  numerical  results  are  converged. 
Otherwise,  any  results  obtained  are  questionable. 

Table  2.5  shows  that  ey  and  e$  only  change  after  the  fourth  decimal  or  at  the 
fifth  significant  figure  with  the  variation  of  C  from  5  to  7.  This  indicates  that  in  the 
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Table  2.5.  Variations  of  «gL,  4JL..  e$»a*,  egL«.  and  Ef]  with  C-  System: 
alternating  trans-polyacetylene;  basis  set:  STO-3G;  geometry:  G2;  energy  unit:  ev. 
For  PLH93,  N=21.  For  PLMBPT,  N=17  and  K=21. 


C=5 


C=6 


C=7 


c(2) 
'c.mm 

.(d) 


e(2) 
M) 

tv.max 


E 
E 


(2) 
9 

(d) 


-0.5779865 
-0.5080903 

0.6403694 
0.5615918 

-1.2183559 
-1.0696818 


-0.5779586 
-0.5080650 

0.6404357 
0.5616481 

-1.2183943 
-1.0697131 


-0.5779526 
-0.5080602 

0.6404428 
0.5616544 

-1.2183954 
-1.0697146 


calculations  of  ey  and  eP' ,  C  =  5  is  already  good  enough  for  current  applications. 
The  variations  of  and  with  C  are  determined  by  those  of  and  e^p  .  Thus, 
the  conclusions  for  E^>  and  E^  are  the  same  as  those  for       and  e"p  . 

Table  2.6  shows  the  variation  of  egL,  e^L.  4JL.,  £j2)  and  with 

Comparing  the  data  for  A'  =  21  and  K  =  41,  one  can  find  that  the  maximum 
difference  for  these  two  groups  is  about  0.4%  as  occurs  for  el2)  .  Thus  for  the 
calculations  where  a  1%  error  is  acceptable,  if  =  21  is  a  reliable  choice.  However, 
when  one  compares  the  data  for  K  -  11  and  those  for  K  —  21,  one  can  find  that  they 
are  quite  different.  The  error  can  be  as  large  as  40%.  This  is  a  very  different  behavior 
from  that  for  E$,  where  only  0.2%  difference  occurred  when  K  was  changed  from 
11  to  21.  Figure  2.1  and  Figure  2.2  show  the  variations  of  Uv>mtx  and  VVtmas  with  ka 
and  kb.  Comparing  to  VVyTnax,  Uv,max  is  flat  and  varies  gradually.  From  Eqs.(2.40)  and 
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Table  2.6.  Variations  of  egL,  e^in,  eg,  „,  e« «,  £f  and  £f  with  K.  System: 
alternating  trans-polyacetylene;  basis  set:  STO-3G;  geometry:  G2;  energy  unit:  ev. 
For  PLH93,  N=21.  For  PLMBPT,  N=17  and  C=7. 


K=ll 

K=21 

K=41 

6(2). 
c,min 

-0.3535667 

-0.5779526 

-0.5809689 

M 

c,min 

-0.3160113 

-0.5080602 

-0.5085070 

e(2) 

0.3954024 

0.6404428 

0.6428650 

M) 

tv,max 

0.3529449 

0.5616544 

0.5614375 

-0.7489691 

-1.2183954 

-1.2238339 

4d) 

-0.6689562 

-1.0697146 

-1.0699445 

(2.50),  we  know  that  E$  only  has  the  terms  like  C/U>max.  Thus,  it  converges  faster 
with  K  while  e„,max  does  not  because  of  Vv,max.  For  e£min,  Vc>min  changes  slowly  with 
ka  and  kb  while  UCiTnin  does  not.  Thus,  e^min  a^so  shows  slow  convergence  with  K. 

Table  2.7  demonstrates  that  c^m,  e{c%in,  eg^ax  and  ei%ax  converge  with  N 
much  slower  than  E$  does.  In  the  calculation  of  E$,  the  difference  between  the 
results  obtained  with  N  =  5  and  N  =  15  is  about  0.57%  while  in  those  for  e^n, 
e*min  ei2mai  and  ei%axi  tne  differences  between  the  results  calculated  with  N=5  and 
N=17  are  larger  than  60%.  The  reason  why  ejj  converges  with  N  so  differently  from 
that  of  E$  is  because  of  the  different  convergence  rate  for  U(P)  and  V(P)  with  N. 
Table  2.8  shows  how  UVtTnax  and  VVtTnax  of  egL*  vary  with  N.  From  the  table,  we  see 
that  Uv<max  converges  very  rapidly  with  N.  The  rapid  convergence  of  UViTnax  with  N 
validate  E^'s  fast  convergence  with  N  since  E$  only  consists  of  a  term  like  UVttnax. 
In  contrast,  Vv>max  shows  slow  convergence  with  N.  Since  Uv>max  and  VViTnax  have 
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Figure  2.2.  Variation  of  the  integrand  in  VViTnax(ka,  kb)  with  ka  and  kb. 
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Table  2.7.  Variations  of  egL,  e{c%n,  e^ax,  eg^,  £f  and  £f  with  N.  System: 
alternating  trans-polyacetylene;  basis  set:  STO-3G;  geometry:  G2;  energy  unit:  ev. 
For  PLH93,  N=21.  For  PLMBPT,  C=7  and  K=21. 


N=5 

N=9 

N=13 

N=17 

eHF 

c.min 
(2) 

c,mm 

3.52581 
-0.94843 

-0.67257 

-0.60361 

-0.57795 

f(d) 
c,min 

-U.oZZyZ 

-u.ooyz  < 

-U.OoUZU 

-U.OUoUO 

,HF  ,A2) 
c,min  '  c,mm 

eHF  +e(d)  ■ 

Z.D  MOO 

2.70289 

Z.o0oZ4 
2.93654 

z.yzzzu 
2.99561 

Z.y4  ( oo 
3.01775 

HF 

v,max 

e(2) 

-4.26684 
1.10012 

0.75796 

0.67279 

0.64044 

0.95061 

0.66210 

0.58944 

0.56165 

eHF  +e(2) 

-3.16673 

-3.50888 

-3.59405 

-3.62640 

eHF  +e(d) 

-3.31623 

-3.60474 

-3.67740 

-3.70519 

7.79265 

5.74411 

6.36212 

6.51625 

6.57426 

6.01912 

6.54128 

6.67301 

6.72294 

opposite  signs  and  e^max  is  the  sum  of  them,  e^max  has  an  even  slower  convergence 
than  K.moi  does. 

2.6    Application  to  Trans  Polvacetvlene 

Trans-polyacetylene  is  probably  one  of  the  most  frequently  investigated  poly- 
mers. One  of  the  topics  of  interest  is  the  behavior  of  the  lowest  transition  energy. 
The  measurement  of  the  absorption  spectrum  of  polyacetylene  crystalline  films  shows 
that  the  absorption  coefficient  begins  a  slow  increase  around  1  eV,  rising  sharply  at 
1.4  eV  to  a  peak  at  about  1.9-2  eV.44  Different  interpretations  exist  for  the  absorption 
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Table  2.8.  Variations  of  UVimax,  K.mai  and  e^max  with  N.  System:  alternating  trans- 
polyacetylene;  basis  set:  ST0-3G;  geometry:  G2;  energy  unit:  ev.  For  PLH93,  N=21. 
For  PLMBPT,  C=7  and  K=21. 


N=5 

N=7 

N=9 

N=ll 

N=13 

N=15 

N=17 

Ui],max 
,(2) 

-1.65723 
2.75734 

1.10012 

-1.65033 
2.51340 
0.86308 

-1.64893 
2.40689 
0.75796 

-1.64881 
2.35252 

0.70371 

-1.64894 
2.32173 
0.67279 

-1.64911 
2.30273 
0.65361 

-1.64933 
2.28978 
0.64044 

spectrum,  namely,  the  assignment  of  the  ~  2  eV  peak  to  the  lowest  singlet  exciton 
band67-72  and  its  assignment  to  a  nir*  interband  transition  where  a  low-energy  tail 
between  ~  1-2  eV  is  then  usually  assigned  to  a  structural  disorder.73-77  A  key  role  in 
this  assignment  problem  is  played  by  the  band  gap  and  the  exciton  binding  energy. 
The  HF  band  gap  is  a  poor  ~6  eV  with  a  large  basis  set.  By  taking  the  second-order 
correlation  into  account  with  MBPT(2),  Suhai  obtained  a  decrease  of  the  theoretical 
band  gap  value  to  about  3  eV.40'41  He  also  showed  that  the  singlet  exciton  band  calcu- 
lated with  correlation  starts  at  about  1.9  eV  and  then  suggested  an  assignment  of  the 
~  2  eV  peak  to  excitons.72  However,  by  solving  the  inverse  Dyson's  equation  in  the 
diagonal  approximation  (correct  to  second-  and  third-order,  respectively),  Liegener 
obtained  a  band  gap  which  is  about  2  eV.39  Then  he  assigned  the  ~  2  eV  peak  to 
the  interband  transition.  Thus,  the  problem  of  how  to  assign  the  ~  2  eV  peak  is  still 
an  unsolved  problem. 

In  this  section,  we  will  compare  our  calculated  total  energies  per  unit  cell  for 
trans-polyacetylene  with  those  obtained  by  Suhai.  Then  we  will  compare  our  results 
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for  the  band  gap  with  Liegener's39  and  Suhai's40'41,  respectively.  Finally,  we  will 
discuss  the  effects  of  the  basis  sets  and  geometries  on  the  band  gap. 

2.6.1    Total  energy  per  unit  cell 

In  this  subsection,  we  will  compare  the  values  of  alternating  trans  polyacetylene's 
E$  calculated  by  PLMBPT  with  those  obtained  by  Suhai40'41.  To  our  knowledge  no 
other  MBPT(2)  data  exist. 

The  numerical  total  energy  per  unit  cell  of  the  system  depends  on  the  geometry 
and  the  basis  set  taken  in  the  calculation.  Table  2.9  lists  the  numerical  values  of  E^CF 
and  E"F  +  obtained  by  Suhai40'41  and  our  results  with  the  same  basis  sets  and 
geometries.  For  the  basis  set  STO-3G  with  the  G2  geometry,  Suhai  did  not  give  the 
value  of  E£)  itself.  We  can  only  compare  the  sum  of  E"f  and  E"F .  They  match 
to  the  first  five  significant  figures,  however.  Considering  that  the  E^F  obtained  by 
PLH93  and  Suhai's  program  agree  in  the  first  five  significant  figures  (see  Appendix 
A),  by  subtracting  the  PLH93's  E%.F  value  from  the  sum  to  estimate  Suhai's  E$. 
In  this  way,  we  estimate  that  the  difference  between  our  results  and  his  is  less  than 
0.4%.  For  the  same  basis  set  STO-3G  but  with  the  G3  geometry,  we  can  directly 
compare  our  MBPT(2)  value  and  his  E$.  From  the  table,  we  see  again  that  the 
difference  between  the  two  results  is  less  than  0.4%.  Considering  there  is  already 
an  appreciable  difference  between  the  band  energies  calculated  by  PLH93  and  those 
given  by  Suhai,  as  indicated  in  Appendix  A,  our  MBPT(2)  results  for  E$  agree  well 
with  Suhai's. 
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Table  2.9.  Comparison  of  E$  calculated  by  PLMBPT  with  those  obtained  by  Suhai. 
System:  alternating  trans-polyacetylene;  energy  unit:  atomic  unit.  For  PLH93, 
N=21.  For  PLMBPT,  K=21,  N=21  and  C=6. 


Basis 

Geo. 

*STO-3G 

G2 

This  worka 

-75.944414 

-0.123390 

-76.067804 

Suhai6 

-76.067437 

STO-3G 

G3 

This  work" 

-75.946823 

-0.118667d 

-76.065490 

Suhaic 

-0.118268d 

-76.065389 

a  E£f  is  calculated  by  PLH93  and  EiV  is  given  by  PLMBPT 
b  Taken  from  reference  41. 
c  Taken  from  reference  40. 
d  Without  two  Is  core  bands 


2.6.2    Comparison  with  Liegener's  results  for  the  band  gap 

As  mentioned  before,  the  band  gap  is  the  difference  between  the  lowest  un- 
occupied band  energy  eC)min  and  highest  occupied  band  energy  ev<max.  For  trans- 
polyacetylene,  eCimin  and  tv>max  appear  at  the  edges  of  the  eighth  and  seventh  bands, 
respectively. 

Table  2  10  rives  the  values  of  e HF     f (d)      e HF  eHF     M)  CHF 

lduie  z.iu  gives  uie  values  oi  ec>min,  ecmin,  ecmin  +  ecmin,  ev  max,  q,^,  ev>max  + 

cKLx.  EgF,  EgF  +  E{gd)  with  a  different  cutoff,  N,  for  the  lattice  summations.  The 

basis  set  is  4-31G**  and  the  geometry  is  G5.  The  second  column  of  the  table  lists 

the  results  obtained  by  Liegener39  with  the  same  basis  set  and  geometry.  They  are 

the  results  obtained  by  solving  the  inverse  Dyson  equation  correct  to  second-order 

in  the  diagonal  approximation  (we  denote  this  method  as  GFM  in  Section  V).  Our 

results  and  his  were  obtained  with  exactly  the  same  formula.  In  his  calculation,  he 
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Table  2.10.  Variations  of  eg!*,  4?L*,  4^«x>  4°  and  4°  with  N  and  com- 

parison  between  our  results  with  Liegener's.  System:  alternating  trans-polyacetylene; 
basis  set:  4-31G**;  geometry:  G5;  energy  unit:  ev.  For  PLH93,  N=21.  For 
PLMBPT,  C=7  and  K=21. 


Liegenera,b 

N=3a 

N=5a 

N=9 

N=13 

N=21 

HF 

c,min 
Ad) 
c,min 

-0.150 
-2.313 

0 

-2 

65579 
68965 

-2.08770 

-1.78252 

-1.69610 

-1.63998 

eHF-  +e(d)  ■ 

-2.463 

-2 

03386 

-1.43191 

-1.12674 

-1.04032 

-0.98419 

v,max 

M) 

-6.400 
2.122 

-5 
3 

78152 
03833 

1.33754 

0.64045 

0.46984 

0.36911 

eHF  +e(d} 

cv ,max  '  cf ,rnax 

-4.278 

-2 

74319 

-4.44398 

-5.14107 

-5.31168 

-5.41241 

6.250 

6 

43731 

1.815 

0 

70933 

3.01207 

4.01533 

4.27136 

4.42822 

a  Without  the  two  Is  core  bands 
b  Taken  from  reference  39 


excluded  the  two  Is  core  bands.  So  for  comparison,  our  results  listed  in  the  third  and 
fourth  columns  were  similarly  obtained.  In  fact,  excluding  the  two  Is  core  bands  in 
MBPT(2)  calculations  only  causes  less  than  a  1%  difference,  which  does  not  affect 
our  current  discussion. 

From  Table  2.10,  one  can  see  that  the  two  HF  band  gaps  are  close  but  there  are 
apprecible  differences  among  the  individual  HF  band  energies,  which  could  influence 
the  second-order  correction  to  the  band  gap.  In  his  paper,  Liegener39  did  not  give  the 
exact  number  of  unit  cells,  N,  included  in  his  calculation.  However,  the  argument  in 
the  paper  suggests  that  the  number  of  the  interacting  cells  was  2,  e.g.  N=3.  From  the 
table,  we  can  also  see  his  results  are  just  between  our  results  obtained  with  N=3  and 


37 


those  with  N=5.  In  Liegener's  calculations,39  Wannier  functions  were  used,  which 
introduce  a  transformation  of  the  wavefunction  and  then  an  extra  transformation  of 
the  two-electron  integrals.  Theoretically,  if  the  transformations  were  done  completely, 
there  would  be  no  difference  between  the  numerical  results  obtained  directly,  and 
those  obtained  via  the  Wannier  functions.  However,  a  large  error  could  occur  if 
the  transformation  were  done  before  convergence  was  reached.  This  could  be  the 
reason  why  his  results  do  not  agree  with  our  results,  but  fell  between  the  two  with 
N=3  and  5,  respectively.  Of  course,  the  number  of  k-points  taken  in  the  calculation 
and  difference  of  the  HF  wavefunction  and  band  structure  could  also  influence  the 
MBPT(2)  results. 

Table  2.10  also  shows  that  the  MBPT(2)  corrections  to  the  band  energies  and 
to  the  band  gap  gradually  approach  converged  results  with  an  increase  in  the  number 
of  unit  cells,  N.  We  already  showed  that  it  is  questionable  to  compare  the  numerical 
results  obtained  with  N  =  3  or  5  with  the  experimental  data,  since  they  are  not 
properly  converged.  Our  final  result  of  the  GFM  band  gap  with  basis  set  4-31G**  and 
geometry  G5  is  4.38270  eV,  contrary  to  Liegener's39  result  of  ~2  eV.  His  assignment, 
based  on  his  calculation,  cannot  be  supported. 

2.6.3    Comparison  With  Suhai's  Results  for  Band  Gap 

Now  comparing  our  results  with  Suhai,  Table  2.11  gives  Suhai's40'41  and  our 
numerical  data  for  four  different  combinations  of  geometry  and  basis  set.  From  the 
table,  one  can  easily  find  large  differences  among  the  two  groups  of  HF  results.  As 
mentioned  in  Section  II,  if  the  cutoffs  in  the  HF  programs  are  well  chosen,  two 
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different  programs  should  give  very  consistent  results.  It  was  shown  in  the  Appendix 
A  that  PLH93  and  CRYSTAL92  provide  very  compatible  results  for  both  the  total 
energy  per  unit  cell  and  the  band  gap  while  their  results,  especially  for  the  band 
gap,  differs  significantly  from  Suhai's.  This  seems  to  suggest  that  Suhai's  HF  results 
were  not  well  converged.  Since  the  HF  results  are  so  different,  it's  hard  to  have  a 
meaningful  comparison  between  our  calculated  MBPT(2)  corrections  to  the  band  gap 
with  his.  It  is  also  not  stated  in  Suhai's  papers  how  many  unit  cells  were  taken,  and 
how  many  k-points  were  chosen  for  the  k-space  integration.  We  know  from  Section 
V.  B  that  the  MBPT(2)  corrections  to  band  energies  are  very  slowly  convergent  with 
both  N  and  K.  Therefore,  this  information  is  very  important  for  judging  the  quality 
of  the  data. 

Our  MBPT(2)  band  gap  with  geometry  G4  and  basis  set  6-31G**  is  3.964  ev, 
which  does  not  support  the  assignment  that  the  measured  ~  2  eV  peak  in  the  trans- 
polyacetylene's  optical  absorption  spectrum  is  the  interband  transition.  We  cannot 
make  such  an  assignment,  before  it  is  certain  that  the  numerical  value  for  the  band 
gap  has  converged  with  basis  set,  and  that  the  geometry  is  accurate,  since  both 
factors  can  have  a  dramatic  effect  on  the  result. 

2.6.4    Effects  of  basis  set  and  geometry  on  the  band  gap 

Table  2.12  shows  how  E«F,  Uv,max,  Vv,max,  e^ax  Uc>min,  Ve>min,  egL  and  E»F+ 
vary  with  four  basis  sets:  STO-3G,  6-31G,  6-31G**  and  DZP.  The  second  column 
is  the  HF  band  gap.  Generally  speaking,  E%F  decreases  with  an  increase  in  the  basis 
set's  size,  but  there  are  exceptions.  For  example,  the  size  of  6-31G  is  smaller  than 
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Table  2.11.  Comparison  of  valence  and  conduction  bands  at  k  =  n/a  and  band  gap  of 
alternating  trans-polyacetylene  between  our  calculated  results  and  Suhai's.  Energy 
unit:  ev.  For  PLH93,  N=21.  For  PLMBPT,  K=21,  N=21  and  C=7.  In  the  MBPT(2) 
calculations,  the  two  Is  core  bands  are  excluded. 


STO-3G  &  G2 

STO-3G  &  G3 

6-31G 

&  G4 

6-31G** 

&  G5 

Suhai"  Ours 

Suhai6  Ours 

Suhai6 

Ours 

Suhai6 

Ours 

eHF 

c,mm 
£(2) 
c,mm 

eHF-  +e(2)  ■ 
eHF 

v,max 

e(2) 

eHF  +e(2) 

EHF 


3.526 

3.719 

3.551 

-0.806 

0.459 

-1.322 

0.576 

-0.563 

-0.340 

-0.584 

-0.430 

-1.274 

-0.674 

-1.902 

3.817 

2.963 

3.379 

2.967 

-1.236 

-0.815 

-1.996 

-1.325 

-4.267 

-4.563 

-4.369 

-5.732 

-5.792 

-5.749 

-5.817 

0.623 

0.270 

0.652 

0.371 

0.749 

0.773 

0.459 

4.604 

-3.644 

-4.293 

-3.717 

-5.361 

-5.044 

-4.976 

-5.358 

7.793 

8.282 

7.921 

4.926 

6.248 

4.427 

6.394 

8.421 

6.607 

7.672 

6.684 

4.125 

4.229 

2.980 

4.033 

°  Taken  from  reference  41 
6  Taken  from  reference  40 


that  of  6-31G**  and  DZP  while  EjfF  calculated  with  6-31G  is  smaller  than  those 
calculated  with  6-31G**  and  DZP.  From  the  third  to  the  eighth  columns,  the  values 

Of  UV!max,  K.maxj  ^v}max  Uc,min,  Ks,min  and  C^min  witn  tne  f°ur  Dasis  sets  are  listed. 

From  Eq.(2.51)  and  Eq.(2.52),  we  know  that  the  number  of  terms  in  both  U(p)  and 
V(p)  increase  with  the  size  of  the  basis  sets  (n).  However,  the  number  of  the  terms  in 
V(p)  is  proportional  to  n  while  that  of  U (p)  is  proportional  to  n2.  It  is  not  necessarily 
true  that  the  more  terms  the  larger  the  magnitude  of  the  sum.  For  example,  VCjTnin 
for  STO-3G  is  the  largest  among  the  four  values.  But  there  could  be  a  tendency 
that  the  more  terms  the  larger  the  magnitude  of  the  sum.  Suppose  this  tendency 
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Table  2.12.  Variations  ofE£F,  Uv,max,  Vv>max  and  eg^  Uc,min,  VCjmin  and  c^in  and 
E^F+E^  with  basis  sets.  System:  alternating  trans-polyacetylene;  geometry:  G6; 
energy  unit:  ev.  For  PLH93,  N=21.  For  PLMBPT,  N=21,  C=7  and  K=21. 


JpHF        JT  T/  A2)  j  j  rr  (2)  J?HFj_T?{1) 

g  uv,max       *v,max    ^v,max        <~» c,min     *c,min        tc,min        g      '  9 


STO-3G  7.12907  -1.67827  2.30916  0.63089  -2.25970  1.70011  -0.55959  5.93859 

6-31G  6.18855  -1.84013  2.58727  0.74713  -2.73089  1.45883  -1.27206  4.16936 

6-31G**  6.35262  -2.43415  2.89380  0.45965  -3.51935  1.61730  -1.90205  3.99092 

DZP  6.28415  -2.45116  3.02154  0.57038  -3.57273  1.60263  -1.97010  3.74367 


exists,  then  the  magnitude  of  U(p)  increases  faster  than  that  of  V(p)  does.  Since  the 
absolute  value  of  Uv>max  is  smaller  than  that  of  the  V^man  the  MBPT(2)  correction 
to  the  highest  occupied  band  energy,  4%ai,  decrease  with  an  increase  in  n  while  for 
the  lowest  unoccupied  band  energy,  the  MBPT(2)  correction,  e^min,  increases  with  n 
because  the  absolute  value  of  f/C)mi„  is  larger  than  that  of  VCtmin.  The  last  column  in 
Table  2.12  lists  the  values  of  the  MBPT(2)  band  gap,  EfF  +  Ef  \  calculated  with 
geometry  G5  and  the  four  basis  sets.  The  numerical  data  show  that  the  larger  the 
basis  set  the  smaller  the  band  gap.  The  smallest  value  among  the  four  is  3.7780  which 
was  obtained  with  basis  set  DZP.  This  value  is  still  much  larger  than  the  ~  2  eV, 
measured  peak  in  the  absorption  spectrum.44  Table  2.12  indicates  that  the  MBPT(2) 
band  gap  could  be  even  smaller  with  larger  basis  sets. 

So  far,  all  the  geometries,  G1-G6,  which  have  been  used  in  our  calculations  were 
optimized  via  the  MBPT(2)  total  energy  per  unit  cell  with  different  basis  sets  by 
Suhai.  The  differences  A  of  the  bond  lengths  rx  and  r2  decrease  with  the  increase 
of  the  size  of  the  basis  set.   The  values  of  A  for  geometries  from  Gl  to  G6  are 
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all  larger  than  0.08 A.  The  geometry  G7  listed  in  Table  2.1  is  a  combination  of 
the  experimentally  estimated  and  theoretical  results.  r3  and  7  are  just  taken  from 
G5.  They  do  not  have  a  significant  effect  on  the  band  gap.  ri,  r2  and  a  were 
experimentally  estimated.66  The  value  of  A  for  G7  is  0.05A,  which  is  much  smaller 
than  all  those  of  the  theoretically  determined  A.  We  know  that  the  theoretical  band 
gap  heavily  depends  on  the  geometry,  especially  A.  When  A  approaches  zero,  we 
obtain  a  symmetric  polyacetylene  which  would  have  a  zero  band  gap.  As  this  point 
corresponds  to  a  different  unit  cell,  and  would  invalidate  non-degenerate  perturbation 
theory,  it  is  outside  the  scope  of  this  study. 

Table  2.13  gives  the  HF  band  gap,  E$ F,  the  MBPT(2)  correction  to  the  band 
gap,  Ef  \  and  the  MBPT(2)  band  gap,  EffF  +  Ef  \  of  the  trans-polyacetylene  with 
five  different  geometries.  The  basis  set  is  DZP  in  all  five  calculations.  The  second 
column  lists  the  results  for  GG  which  has  A  =  0.0816.4.  The  third  column  lists 
the  results  for  G7,  the  A  of  which  is  0.05A.  Starting  from  the  fourth  column,  the 
structural  parameters  r3,  a,  7  and  ri+r2  are  fixed  to  be  the  same  as  those  for  G7,  and 
only  A  changes.  The  results  in  Table  2.13  show  that  the  HF  band  gap  EffF  decreases 
with  a  decrease  in  A,  as  does  the  MBPT(2)  correction.  The  MBPT(2)  band  gap 
reduces  to  2.1884  eV  when  A  is  decreased  to  0.004A.  Of  course,  the  MBPT(2)  band 
gap  depends  on  other  structural  parameters  such  as  the  length  of  the  unit  cell  and 
q.  Since  there  is  a  big  difference  between  the  theoretically  determined  geometry  and 
the  experimentally  estimated  geometry  upon  which  the  band  gap  heavily  depends, 
further  investigation  on  the  geometry  of  trans-polyacetylene  with  larger  basis  sets  are 
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Table  2.13.  Variations  of  E%F,EW  and  E^F+E^  with  A  =  r2  -  rx.  System: 
alternating  trans-polyacetylene;  Basis  set:  DZP;  Energy  unit:  ev.  For  PLH93,  N=21. 
For  PLMBPT,  N=21,  C=7  and  K=21. 


Geometry 

G6 

G7 

A  =  0.03A 

A  =  O.OlA 

A  =  0.004A 

E{g2) 

6.28415 
2.54048 
3.74367 

5.57441 
2.35208 
3.22233 

5.08259 
2.28585 
2.79674 

4.53567 
2.20086 
2.33481 

4.35645 
2.16935 
2.18710 

needed.  Since  the  band  gap  is  a  sensitive  quantity,  higher  order  corrections  may  have 
significant  contributions  too.  The  difference  between  the  real  experimental  condition 
(crystallized  film  of  trans-polyacetylene  and  even  differences  due  to  different  ways 
of  processing  the  polymers)  and  the  theoretical  model  of  an  isolated  infinite  chain 
may  be  another  reason  for  the  mismatch  between  the  MBPT(2)  band  gap  and  the 
measured  ~2  eV  peak  in  the  absorption  spectrum.44 

2.7  Conclusions 

In  this  chapter,  we  have  presented  explicit  expressions  of  MBPT(2)  for  the 
total  energy  per  unit  cell  and  explicit  expressions  for  a  direct  evaluation  of  the  band 
structure  and  implemented  them  for  infinite,  periodic  one  dimensional  polymers. 
Three  independent  programs  are  used  to  establish  the  veracity  of  the  results.78  The 
convergence  behavior  of  MBPT(2)  corrections  both  to  total  energy  per  unit  cell  and  to 
band  energies  have  been  studied.  It  is  found  that  the  MBPT(2)  correction  converges 
very  fast  with  the  number  of  unit  cells  and  the  number  of  k-points  taken  for  the  k- 
space  integration  over  the  first  Brillouin  zone,  while  the  MBPT(2)  correction  to  the 
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band  structure  converges  much  slower  with  N  and  requires  a  much  larger  K.  Neither 
MBPT(2)  corrections  are  sensitive  to  an  integral  cutoff  of  10c  when  C  >  5. 

In  application  to  alternating  trans-polyacetylene,  the  DZP-MBPT(2)  numerical 
values  of  the  band  gap  are  larger  than  the  measured  ~2  eV  peak  in  the  absorption 
spectrum  of  the  system.  Further  investigation  with  larger  basis  sets  and  better  op- 
timized geometries,  and  higher-order  correlation  corrections  are  required  to  make  an 
assignment  of  the  measured  ~2  eV  peak  either  as  an  exciton  band  or  as  an  inter- 
band  transition.  Future  work  will  address  higher-order  correlation  corrections  with 
coupled-cluster  theory  and  application  to  systems  periodic  in  two  and  three  dimen- 
sions. 


CHAPTER  3 

THE  MBPT(2)  PHOTOELECTRON  SPECTRUM  OF  POLYETHYLENE: 
EXPLANATION  OF  XPS  AND  UPS  MEASUREMENTS 


Polyethylene  is  probably  the  best  studied  polymer  by  XPS  and  UPS.  The  mea- 
sured XPS  and  UPS  spectra  provide  rich  information  about  the  electronic  structure 
of  the  system  and  a  sensitive  to  tool  to  check  theory. 

Various  semiempirical  methods52  and  Hartree-Fock  (HF)  calculations79-81  have 
been  used  to  attempt  to  intepret  polyethylene's  UPS  and  XPS  spectra.  A  few  semiem- 
pirical methods,  among  many,  provided  some  agreement  with  experiment  in  selected 
energy  ranges,  but  failed  in  others.52  HF  offered  a  better  description  of  the  general 
features,79-81  but,  in  the  absence  of  electron  correlation,  which  is  crucial  in  a  first 
principle  description  for  both  finite  molecules  and  extended  systems,1'5'10-12  the  HF 
band  structure  would  require  a  2  eV  shift  and  an  80%  contraction  in  scale  to  compare 
with  experiment.51  However,  there  are  discrepancies  in  the  experiments  that  requires 
resolution.  The  experimental  XPS  data47  finds  a  peak  at  12.6  eV  that  is  not  seen  in 
either  the  other  XPS  work48  or  ARUPS  experiments,52  and  misses  a  peak  at  15.4  eV, 
seen  in  the  other  experiments.48-52  The  ARUPS  result  misses  the  9.6  ev  XPS  peak.52 
Predictive  correlated  theory  can  help  to  resolve  these  questions. 


44 


45 


In  this  chapter,  we  report  the  initial  MBPT(2)  quasi-particle  valence  band  ener- 
gies for  polyethylene.  We  will  show  that  unlike  DFT,  the  ab  initio  MBPT(2)  descrip- 
tion in  a  polarized  6-31G**  basis,82  properly  converged  with  lattice  summations  and 
integration  over  reciprocal  space,43  accurately  explains  the  XPS  and  UPS  spectra. 
We  will  also  show  that  the  disagreements  among  the  experiments  can  be  resolved  by 
considering  the  width  of  the  X-ray  and  UV  radiation. 

Polyethylene  is  considered  in  its  all-  trans  conformation  which  has  a  screw  axis  52 
along  the  chain  direction.  The  unit  cell  is  CH2.  Unlike  polyacetylene,  there  are  only 
slight  differences  among  the  geometries  determined  either  by  experiment  or  by  theory. 
In  our  calculation,  we  will  use  the  X-ray  structure:  rCc  =  2.89  bohr,  rCn  —  2.02 
bohr,  IHCH  =  107.0  deg  and  ICCC  =  112  deg.80-83 

Our  MBPT(2)  program  for  extended  systems  for  the  total  energy  and  band 
structure  is  combined  with  HF  solutions  from  PLH93.30'43  In  our  calculations,  we 
employ  34  unit  cells  (CH2)  in  the  lattice  summations  and  41  points  in  the  first  Bril- 
louin  zone  for  integration  over  the  reciprocal  space  to  ensure  that  the  numerical  results 
are  converged.  This  is  critical  for  MBPT(2)  calculations  in  extended  systems.43 

Figure  3.1  shows  that  the  HF,  two  DFT  variants  (LDA  and  gradient  corrected 
BLYP)84,  and  MBPT(2)  band  structures  of  polyethylene  calculated  with  the  6-31G 
basis  set.  The  MBPT(2)  bands  are  above  the  HF  and  below  the  DFT  ones.  The 
correlation  shift  is  different  at  different  points  in  the  bands,  being  around  2  eV  for 
the  first  two  bands.  For  the  third  band,  the  shift  is  about  5  eV  at  0  and  3  eV  at  n/a. 
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Figure  3.1.  The  HF,  LDA,  BLYP,  and  and  MBPT(2)  valence  bands  of  all-trans 
polyethylene  with  basis  set  6-3 1G  as  a  function  of  k. 

Figure  3.2a  shows  the  density  of  states  (DOS)  for  the  MBPT(2)  bands  shown  in 
Figure  3.1.  The  peaks  in  the  DOS  have  corresponding  maxima  in  the  photoelectron 
spectra  and  can  be  well  determined  experimentally.  The  calculated  values  for  these 
peaks  using  HF  and  MBPT(2)  with  6-31 G  and  6-31G**  bases,  that  latter  containing 
d  functions  on  C  and  p  functions  on  H,  are  compared  to  experiment  in  Table  1.  As 
polarization  functions  are  essential  to  measuring  correlation  effects,  we  see  shifts  of 
0.1  to  0.3  eV. 

The  XPS  and  UPS  for  polyethylene  are  measured  in  the  solid  phase.  However, 
it  is  known  that  the  weak  van  der  Waals  intermolecular  interaction  causes  only  about 
a  0.1  eV  energy-band  dispersion,85  which  is  negligible  compared  with  the  large  in- 
tramolecular dispersion  of  5  eV.  The  difference  between  the  binding  energy  of  an 
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Figure  3.2.  The  MBPT(2)  photoelectron  spectra  for  all-trans  polyethylene  with  basis 
set  6-31G.  (a)  density  of  states;  (b)  photoelectron  spectrum  with  F  =  0.2  eV;  (c) 
photoelectron  spectrum  with  F  =  0.75  eV. 
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Table  3.1.  Comparison  among  the  peaks  in  the  density  of  states  calculated  using  HF, 
BLYP,  and  MBPT(2)  with  basis  sets  6-31G  and  6-31G**,  respectively,  and  those 
measured  by  XPS  and  ARUPS. 


LP. 

Ii 

I2 

I, 

I3 

I4 

Hi 

II2 

HF 

6-31G 

10.5 

12.31 

12.52 

16.08 

18.37 

21.43 

30.21 

HF 

6-31G** 

10.6 

12.37 

12.64 

15.99 

18.27 

21.50 

30.03 

BLYP 

6-31G 

5.9 

7.39 

7.40 

10.24 

11.73 

13.70 

19.92 

MB6 

6-31G 

8.1 

9.5 

10.32 

10.75 

13.44 

15.91 

18.48 

25.09 

MB6 

6-31G** 

8.4 

10.59 

11.06 

13.54 

15.99 

18.39 

24.66 

XPS47 

PE6 

8.6 

9.6 

11.2 

12.6 

13.8 

18.0 

23.6 

XPS48 

C36H74 

9.8 

11.1 

13.7 

15.4 

18.0 

23.8 

UPS8 

C36H74 

10.5-12.0C 

14.0 

15.5 

18.3 

24.6 

a  MBPT(2)/6-3lG**  calculations  have  been  done  only  at  the  peaks. 
b  MB  and  PE  denote  MBPT(2)  and  polyethylene,  respectively. 
c  depends  on  the  parameters  of  ARUPS  measurements.52 


individual  chain  in  the  gas  phase  and  that  in  the  solid  phase  is  the  work  function, 
a  constant  determined  to  be  in  the  range  of  4.5  to  4.8  eV.  The  experimental  values 
listed  in  Table  1  take  the  work  function  as  4.8  eV.52 

The  line  intensities  of  the  measured  spectra  depend  on  both  the  photoionization 
cross  sections  and  the  frequency  distribution  of  the  incident  radiation.  The  photoion- 
ization cross  sections  are  functions  of  the  angle  of  incident  radiation,  the  energy  of  the 
radiation  and  the  angle  of  the  emitted  electrons.  Hence,  the  ARUPS  spectra  varies 
with  these  three  parameters.52  Here,  we  only  focus  on  the  peaks  of  the  DOS  when  we 
compare  with  the  ARUPS  spectra.  For  the  XPS  spectra,  we  simply  use  the  Gelius 
model86  for  the  relative  photoionization  cross  section  taking  the  photoionization  cross 
section  of  the  C2s  bonding  band  to  be  13  times  larger  than  that  for  the  other  two 
valence  bands. 
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The  energy  distribution  of  the  incident  radiation  can  be  described  by  a  linear 
combination  of  Lorentzian  and  Gaussian  curves  with  a  half  width  at  half-maximum 
T.  Then  the  relative  line  intensities  are  expressed  as 


dE',(3.1) 


where  tt  is  the  ith  band's  photoelectron  cross  section,  di(E')  is  the  band's  DOS,  while 
wi  and  wg  are  the  weight  of  the  Lorentzian  and  Gaussian  curves,  respectively. 

With  T  =  0.2  eV  and  using  equal  weight  for  the  two  curves,  we  obtain  the 
MBPT(2)  photoelectron  spectrum  described  in  Figure  3.2b.  Because  of  the  strong 
overlap,  L  and  I2  merge  into  one  peak  as  do  the  two  peaks  of  Hi .  Other  peaks  in 
the  DOS  have  their  own  corresponding  maxima  in  Figure  3.2b.  The  first  ionization 
potential  (LP.)  has  a  small  but  visible  peak. 

The  full  width  at  half-maximum  of  the  X-ray  used  in  the  XPS  experiments  was 
about  1.5  eV,87  e.g.  T  =  0.75  eV.  Using  Eq.(3.1)  with  equal  weight  for  the  two 
curves,  we  obtain  the  corresponding  line  intensities  described  in  Figure  3.2c.  There 
is  no  other  peak  corresponding  to  the  first  ionization  potential  in  Figure  3.2c  since 
its  photoelectron  cross  section  is  diffuse.  This  agrees  well  with  both  the  XPS  and 
UPS  spectra,  in  which  no  peak  around  the  LP.  in  Figure  3.2a  was  observed. 

The  measured  shoulder  in  XPS,  labeled  Is  in  Table  1,  was  originally  erroneously 
thought47-80  to  be  the  LP.,  but,  instead,  is  the  result  of  overlap  among  the  Ilf  I2  and 
the  LP.  peaks.  The  enlarged  part  in  Figure  3.2c  clearly  shows  that  there  is  a  shoulder 
in  the  MBPT(2)  spectrum  when  L  is  0.75  eV.  The  shoulder  located  at  about  9.5  eV 
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agrees  well  with  the  XPS  results  of  9.6  eV  and  9.8  eV.  The  shoulder  appears  in  the 
photoelectron  spectrum  only  when  the  relative  photoelectron  cross  sections  of  the  two 
bands  corresponding  to  C2p  and  HXs  bonding  have  a  suitable  ratio.  Because  of  this 
unfavorable  ratio,  the  shoulder  could  not  be  seen  in  most  UPS  spectra.52  However,  if 
one  checks  the  measured  spectra  carefully,  one  can  discern  the  shoulder  in  the  spectra 
observed  at  special  angles  and  energies. 

Althought  the  LP.  cannot  be  accurately  measured  in  either  the  XPS  or  UPS 
experiments,  since  the  radiation  width  is  not  small  enough,  it  can  be  deduced  from 
the  measured  energy  gap  and  the  work  function.  The  LP.  value  determined  by 
Delhalle  et  al47  in  this  way  for  polyethylene  was  8.3  eV.  The  work  function  used  in  his 
determination  was  4.5  eV.  To  match  with  the  other  experimental  data  given  by  Seki 
et  al.,52  the  experimental  IP,  namely  8.6  eV,  listed  in  Table  1  is  determined  by  taking 
the  work  function  to  be  4.8  eV  deduced  from  the  ARUPS  experiments.  Our  MBPT(2) 
I.P.with  basis  set  6-31G**  is  8.4  eV,  which  agrees  well  with  the  experimental  result. 

In  both  the  XPS  and  UPS  spectra,  the  two  peaks  corresponding  to  Ii  and  I2  are 
too  close  to  be  distinguished,  so  only  one  peak  is  observed  with  any  given  experimental 
parameters.  As  mentioned  for  I,  above,  the  relative  photoelectron  cross  sections  of 
the  two  peaks  vary  differently  with  the  incident  angle  of  the  radiation,  the  energy  of 
the  radiation,  and  the  emission  angle  of  the  electrons.  The  maximum  corresponding 
to  Ix  and  I2  in  the  photoelectron  spectra  varies  between  Ix  and  I2  with  the  three 
parameters.  It  could  be  larger  than  I2  since  the  DOS  at  the  outside  of  I2  is  very 
large.  The  position  of  the  peak  measured  in  the  ARUPS  experiment,  indeed,  varies 
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from  10.5  to  12.0  eV  with  the  experimental  parameters.52  Considering  our  MBPT(2) 
values  for  Ix  and  I2  with  basis  set  6-31G**  are  10.59  and  11.06  eV,  respectively,  the 
agreement  between  theory  and  experiment  is  excellent.  Since  the  three  parameters  in 
the  XPS  experiment  are  fixed,  only  one  peak  is  observed  for  Ii  and  I2,  which  occurs 
at  11.1  and  11.2  eV  in  the  two  XPS  measurements,  respectively,  again  in  excellent 
agreement. 

The  peak  1^  measured  in  an  earlier  XPS  experiment47  was  not  observed  in 
the  later  XPS  spectra48  and  the  ARUPS  measurements.51  The  MBPT(2)  ionization 
spectrum  in  Figure  3.2c  supports  the  later  experiments,  indicating  the  lx  peak  to  be 
an  artifact. 

For  the  peaks  I3  and  I4,  the  MBPT(2)  values  agree  closely  with  those  observed 
in  the  three  experiments.  The  peak  I4,  as  shown  in  Figure  3.2c,  is  at  the  foot  of  the 
much  stronger  peak  Hi.  It  is  easily  suppressed  by  the  tail  of  Hi  if  the  signal-to-noise 
ratio  is  not  adequate.  This  is  probably  why  I4  was  not  observed  in  the  earlier  XPS 
experiment47  while  it  was  in  others.48'51 

The  two  energetic  peaks  Hi  and  II2  come  from  the  C2s  bonding  band.  Hi  has 
two  peaks  in  the  spectrum  of  the  DOS  described  in  Figure  3.2a.  The  peaks  are 
so  close  that  only  one  peak  could  be  observed  in  the  XPS  and  UPS  experiments, 
as  happened  for  the  two  peaks  corresponding  to  Ix  and  I2.  However,  since  the  two 
peaks  of  Hi  belong  to  the  same  band,  the  relative  photoelectron  cross  sections  of  the 
two  peaks  vary  in  the  same  pattern  with  the  two  angles  and  the  radiation  energy 
in  the  ARUPS  measurements.  Thus,  the  position  of  the  peak  should  not  vary  much 
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with  these  experimental  parameters.  This  is  why  a  peak  with  a  stable  position  for 
Hi  was  observed  while  the  peak  for  Ii  and  I2  varied  over  a  range  in  the  ARUPS 
measurements.52  Once  again,  our  MBPT(2)  results  agree  well  with  those  observed  in 
the  XPS  and  ARUPS  experiments. 

For  the  II2  peak,  the  three  experimental  values  are  23.6,  23.8  and  24.6  eV, 
respectively,  differing  by  1  eV.  Our  MBPT(2)/6-31G**  value  is  24.66  eV,  which  is  in 
excellent  agreement  with  that  of  the  ARUPS  measurement,  believed  to  be  the  most 
reliable.52 

A  direct  comparison  between  the  MBPT(2)  photoelectron  spectra  shown  in 
Figure  3.2c  and  the  experimental  XPS  spectra  measured  by  J.  J.  Pireaux,  et  al486  is 
given  in  Figure  3.3,  where  the  Hi  peak  of  the  two  spectra  have  been  superimposed. 
The  two  spectra  match  very  well,  even  including  the  three  small  peaks  that  fall 
between  -10  and  -16  on  Figure  3.3. 

In  Figure  3.4,  we  compare  the  simulated  MBPT(2)  spectrum  with  that  from  HF 
and  DFT  (LDA  and  BLYP  are  indistinguishable).  Clearly,  only  MBPT(2)  properly 
describes  the  spectrum.  HF  shows  the  three  small  peaks,  but  substantially  shifted, 
while  DFT  seems  to  show  two  of  the  features,  the  third  being  lost  in  the  first  large 
peak.  No  constant  shift  will  match  the  HF  or  DFT  spectra  to  MBPT(2). 

In  conclusion,  we  have  obtained  the  quasi-particle  band  energies  for  the  three 
valence  bands  of  polyethylene  using  correlated  (MBPT(2),  6-31G,  6-31G**)  ab  initio 
theory  for  extended  systems.  Unlike  HF  or  DFT  (LDA  and  BLYP)  calculations,  the 
MBPT(2)  explain  the  XPS  and  UPS  measurements  accurately,  while  also  resolving 


53 


-30  -20  -10  o 

Figure  3.3.  Comparison  between  the  experimental  and  the  MBPT(2)  XPS  for 
polyethylene.  Solid  line:  experimental  XPS  measured  by  J.  J.  Pireaux  et  al486;  Cir- 
cles: the  MBPT(2)  photoelectron  spectra  shown  in  Fig.  2c. 
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Figure  3.4.  Comparison  between  the  MBPT(2)  spectra,  with  that  from  HF  and  DFT. 
The  DFT  results  from  either  LDA  and  gradient  corrected  BLYP  are  indistinguishable 
on  the  figure. 

the  discrepancies  among  the  experiments.  Further  work  will  focus  on  higher  order 
CC/MBPT  and  two  and  three  dimensional  systems. 
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CHAPTER  4 

MBPT(2)  CALCULATIONS  FOR  POLYMETHINEIMINE:  GEOMETRY,  BAND 
GAP,  AND  VIBRATIONAL  FREQUENCIES 


4.1  Introduction 

Polymethineimine  was  first  synthesized  in  the  early  seventies  by  Wohrle.54  It  is 
the  prototype  of  a  class  of  polymers  which  possess  conjugated  carbon-nitrogen  double 
bonds.  Conductivity  measurements  indicate  that  polymethineimine  is  a  semiconductor55, 
making  it  a  potential  organic  conductor  with  doping.88  It  also  has  interesting  non- 
linear optical  properties  because  it  has  an  asymmetric  monomer.89  However,  except 
for  its  three  vibrational  frequencies,54'55  little  about  this  polymer  has  been  measured. 

Karpfen90  and  Teramae  et  al91  have  calculated  the  stable  structure  for  polyme- 
thineimine using  the  Hartree-Fock  (HF)  method  with  basis  sets  6-31G  and  STO-3G, 
respectively.  Their  results  confirmed  the  conjugated  structure  of  the  polymer,  but 
there  were  obvious  numerical  differences  among  the  two  calculations  since  the  basis 
sets  used  in  their  calculations  were  not  large  enough  to  obtain  converged  results,  nor 
did  their  calculations  include  electron  correlation  which  is  likely  to  be  critical  in  de- 
termining the  equilibrium  structures  and  other  properties,  especially  for  unsaturated 
systems. 
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For  the  band  gap,  empirical,92  semi-empirical,93  and  HF91  methods  have  been 
applied  to  polymethineimine,  with  HF  theory  giving  the  largest  band  gap.  From  the 
HF  calculation  mentioned  above,  Teramae  et  al91  conclude  that  polymethineimine 
was  unlikely  to  be  an  organic  conductor  with  doping.  This  contradicts  the  observed 
semiconducting  behavior  of  the  system.91 

Vibrational  frequencies94  of  a  system  provide  important  information  about  the 
system's  structure,  as  they  are  very  sensitive  to  it.  They  become  even  more  impor- 
tant when  the  equilibrium  structure  is  not  known.  Vibrational  frequencies  are  also 
important  properties  for  either  finite  molecules94  or  infinite  chains.95  Since  they  are 
determined  by  the  second  derivatives  of  potential  energy  surfaces,  they  are  sensitive 
to  the  accuracy  of  the  theoretical  method  and  the  size  of  basis  set.12  Few  vibrational 
frequency  calculations  using  ab  initio  methods  have  been  done  for  infinite  chains,5,91 
and  all  of  those  calculations  were  performed  at  the  HF  level  or  lower.5'91 

Teramae  et  al91  obtained  the  frequencies  for  polymethineimine  using  HF  in  a 
STO-3G  basis  set,  scaling  with  a  factor  of  0.89.  However,  lacking  electron  correlation 
and  basis  set  flexibility  in  their  calculation,  little  can  be  said  about  the  convergence 
of  the  calculation. 

In  this  chapter,  we  apply  MBPT(2)  to  polymethineimine  to  determine  its  equi- 
librium structure,  its  band  gap,  and,  particularly,  its  vibrational  frequencies.  Our 
objection  is  to  assess  the  effect  of  the  simplest  level  of  correlation  on  the  phonon 
spectra,  since  no  such  prior  calculation  has  been  made. 
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The  plan  of  the  paper  is  as  follows.  In  section  4.2,  we  will  briefly  describe  the 
procedure  for  the  geometry  optimization  in  extended  systems,  then  we  will  give  the 
basic  formulas  for  the  band  gap  calculation,  and  finally  we  will  discuss  the  calculation 
of  vibrational  frequencies  in  infinite,  periodic  systems.  In  section  4.3,  we  will  apply 
the  theory  described  in  Section  4.2  to  polymethineimine  to  determine  its  equilib- 
rium structure,  its  band  gap,  and  vibrational  frequencies.  Then  we  will  follow  with 
conclusions. 

4.2  Theory 

4.2.1    Geometry  optimization 

For  a  one-dimensional  periodic  system,  the  independent  variables  to  describe  the 
nuclear  structure  are  the  internal  coordinates  in  one  unit  cell  and  the  unit  cell  length 
a.  If  the  symmetry  of  the  system  is  larger  than  just  the  translational  symmetry,  the 
number  of  independent  coordinates  can  be  further  reduced.  Using  r  =  n,  r2, rm  to 
denote  the  independent  internal  coordinates  of  the  polymer,  then  the  Hartree-Fock 
(HF)  total  energy  per  unit  cell  Ej£ f)  and  the  MBPT(2)  total  energy  per  unit  cell 

EMBPT(2)  for 

a  given  basis  set  are  functions  of  rx,  r2,  rm. 
embpt(2)  ig  the  sum  of  E{hf)  and  the  Second-order  MBPT  correction  to  the 
total  energy  per  unit  cell,  E$,  e.g. 

Eg*'™  =  Eg},  (4.1) 

where  EgF  and  E$  can  be  calculated  by  Eqs.  (2.19)  and  (2.40),  respectively. 
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The  optimized  HF  and  MBPT(2)  geometries  with  the  basis  set  are  the  geome- 
tries at  which  E^F^  and  Ef?cBPTW  have  minimum  values,  respectively.  There  are 
many  efficient  algorithms  developed  to  locate  minima  on  potential  energy  surfaces.96'97 
Among  them,  the  quadratic  steepest  descent  line  following  algorithm98  is  more  ef- 
ficient and  stable.  The  basic  idea  of  this  method  is  to  construct  a  quadratic  form 
using  the  energy  and  the  gradient  at  the  current  point  with  an  updated  hessian,  then 
follow  the  steepest  descent  line  of  the  quadratic  form  to  the  next  point,  and  repeat 
this  procedure  until  convergence  is  reached. 

Optimized  geometries  or  minima  on  potential  energy  surfaces  correspond  to 
molecular  equilibrium  structures.  Besides  their  importance,  they  also  provide  the 
starting  point  for  the  determination  of  all  other  properties. 

4.2.2    Band  gap 

In  an  extended  system,  orbital  energies  enk  form  the  energy  bands.  According 
to  Koopmans'  theorem,63  the  HF  band  energies  for  qusi-particle  occupied  and  unoc- 
cupied bands  are  equal  to  the  corresponding  HF  ionization  potentials  and  electron 
affinities,  respectively.  Beyond  the  HF  approximation,  one  can  still  define  the  band 
energies  as  the  ionization  potentials  for  occupied  orbitals  and  the  electron  affinities 
for  unoccupied  orbitals.37 

The  MBPT(2)  orbital  energies  are  defined  as  the  MBPT(2)  energy  difference 
between  the  neutral  and  the  ion  systems  and  are  given  by 


MBPT(2)  _    (HF)  (2) 

tP  —  eP     +eP  , 


(4.2) 
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where  e(P'  is  the  second-order  MBPT  correction  to  the  band  energy  and  can  be 
calculated  by  using  Eqs.  (2.47),  (2.50),  and  (2.51). 

The  band  gap  in  an  extended  system  is  defined  as  the  difference  between  the 
lowest  band  energy  on  the  conduction  band  and  highest  band  energy  on  the  valence 
band.37  It  determines  whether  the  system  is  a  conductor,  a  semiconductor  or  an 
insulator.  The  HF  band  gaps  are  much  larger  than  the  measured  values  due  to 
the  lack  of  electron  correlation.  Including  electron  correlation,  MBPT(2)  band  gaps 
greatly  improve  the  HF  results  although  they  are  still  larger  than  the  experimental 
values,  primarily  because  of  the  corresponding  poor  description  of  the  conduction 
bands. 

4.2.3    Vibrational  frequencies 

Piseri  and  Zerbi"  have  developed  the  theory  for  vibrational  motions  of  nuclei 
in  an  infinite  periodic  chain."  In  the  following,  we  will  give  the  basic  formulas  to  be 
used  in  our  calculations. 

Let  r„  and  xn  be  the  displacements  of  all  internal  and  Cartesian  coordinates, 
respectively,  relative  to  the  stable  structure  in  the  nth  unit  cell.  Hence,  we  have  the 
relationship 


rn  —  XI  B/Xn+/, 


(4.3) 


60 


where  Bj  are  transformation  matrices.  Then  the  kinetic  energy  and  potential  energy 
for  the  system  can  be  expressed  as 


r=^xtMxn,  (4.4) 


where  M  is  the  mass  matrix  which  is  diagonal,  and 


V'=^ErnF""'rn''  (4-5) 


where  Fnn>  are  the  force  constants  with  respect  to  rn  and  r„/,  respectively.  The 
translational  symmetry  ensures  that 

Fnn'  =  F0n'_„.  (4.6) 

Using  R(fc)  and  X(/c)  to  denote  the  symmetrized  internal  and  Cartesian  coordinates 
for  translational  symmetry.  Then  TL(k)  and  X(k)  can  be  calculated  by 


R(A;)  =  ^==^rne-fc"/a  (4.7) 
J2n/a  n 


and 


X(k)  =  -rL=J2xne-^a,  (4.8) 
J2ir/a  n 
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respectively.  The  internal  and  Cartesian  coordinates  for  the  symmetried  coordinates 
~R(k)  and  X(A;)  can  also  be  calculated  by 


rn  =         —^eikn'aR(k)dk,  (4.9) 


xn  =         —L=eikn/aX(k)dk.  (4.10) 

J-ir/a  .  /9-7T  In 


Substituting  Eqs.  (4.9)  and  (4.10)  into  Eqs.  (4.4)  and  4.5),  one  can  get 


T=  \        X{k^MX(k)dk,  (4.11) 

2  J—ir/a 


V=\  [W/a  R(kyFR(k)R(k)dk,  (4.12) 

2.  J—nja 


where 


FR(k)  =  Y,Fonelkn/a  (4.13) 


and  use  has  been  made  of  Eq.(4.6)  and 


n 


(4.14) 
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From  Eqs.  (4.3),  (4.9),  and  (4.10),  we  can  connect  R(/c)  and  X(A;)  by 

R(fc)  =  B(ifc)X(fc)  (4.15) 

where 

B(k)  =  J^e-^Bi.  (4.16) 
Substituting  Eq.(4.15)  into  Eq.(4.12),  we  get 

V  =  \  f'a  X(k)^Fx(k)X(k)dk,  (4.17) 

L  J —it /a 

where 

Fx  (*)  =  B(kyFR(k)B(k).  (4.18) 

Then  the  vibrational  frequencies  co(k)  can  be  calculated  by 

Det\M~1/2Fx(k)M~1/2  -  u{k)I\  =  0,  (4.19) 

where  I  is  an  identity  matrix.  The  rank  of  Fx(k)  is  determined  by  the  number  of 
Cartesian  coordinates  in  one  unit  cell,  providing  the  same  number  of  solutions  for  a;  (A:) 
from  Eq.(4.19).  If  there  is  an  additional  symmetry  and  only  the  vibrational  modes  for 
certain  representations  are  of  interest,  the  size  of  the  matrix  Fx(k)  becomes  smaller. 
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The  solutions  of  Eq.(4.19)  form  the  phonon  disperson  curves  which  are  analogous  to 
energy  bands.  From  the  phonon  disperson  curves,  the  density  of  the  phonon  states 
can  be  determined.  From  Eqs.  (4.9)  and  (4.10),  we  can  see  that  for  a  nucleus  moving 
with  a  phase  <j>,  there  is  always  a  corresponding  nucleus  in  another  unit  cell  moving 
with  a  phase  of  —0  in  a  vibrational  mode  with  k  ^  0  or  n/a.  Their  contributions 
to  the  dipole  momentum  always  cancel  each  other.  Then  all  vibrational  motions  for 
k  ^  0  or  ir/a  do  not  change  the  dipole  moment.  The  corresponding  vibrational  modes 
are,  therefore,  not  active  and  cannot  be  observed  in  the  IR  spectra.95  For  k  =  0,  all 
the  corresponding  nuclei  in  different  unit  cells  move  in  phase  and  their  contributions 
to  the  dipole  moment  add.  Since  the  value  of  the  contribution  from  each  unit  cell  is 
a  function  of  the  nuclei'  coordinates,  the  dipole  moment  changes  with  the  vibrational 
motions.  Then  the  vibrational  modes  with  k  =  0  are  observable  in  an  IR  spectra. 
The  frequencies  for  these  modes  are  called  fundamental  frequencies. 

From  Eq.(4.9),  the  translational  symmetry  is  preserved  when  k  =  0  and  then 

ehf  ,  E{2)    _    I  VrtF.  r 

^  71 

=   2rtPr>  (4-20) 

where  r  are  internal  coordinates  in  any  one  unit  cell.  Taking  the  second  derivatives 
on  both  sides  of  Eq.(4.20),  we  get 

<F*<0»o'  =  |:J:[££F+^]-  (4.2i) 

1  3 
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For  a  polymer  with  screw  symmetry  such  as  polyethylene,  the  vibrational  modes 
for  k  =  n/a  may  also  be  observed.  Otherwise,  the  vibrational  modes  for  k  —  ir/a  are 
not  active  in  IR  spectra. 

4.3    Application  to  Polymethineimine 

PLH9332  and  PLMBPT43  programs  are  used  for  HF  and  MBPT(2)  calcula- 
tions, respectively.  The  number  of  unit  cells  for  lattice  summations  in  both  HF  and 
MBPT(2)  calculations  is  21.  Multipole  expansions  have  been  used  in  HF  calculations 
for  the  lattice  summations  outside  of  the  21  unit  cells  to  further  enhance  the  con- 
vergence of  the  HF  results.32  The  convergence  of  MBPT(2)  for  both  the  total  energy 
per  unit  cell  and  for  the  quasi-particle  band  energies  with  the  lattice  summations 
have  been  studied  in  our  previous  report.43  The  total  energy  per  unit  cell  converges 
much  faster  with  the  lattice  sums  than  the  quasi-particle  band  energies.  For  our  cur- 
rent application,  21  unit  cells  in  the  lattice  summations  are  sufficient  to  give  reliable 
numerical  results. 

4.3.1    Stable  structure 

The  structure  of  polymethineimine  has  not  been  determined  by  experiment  yet. 
There  are  a  few  theoretical  geometries  for  this  polymer  obtained  by  HF  calculations 
with  small  and  medium  sized  basis  sets.90,91  In  the  following,  we  will  determine  the 
geometry  using  both  HF  and  MBPT(2)  with  different  basis  sets. 

Figure  4.1  shows  the  structure  of  all-trans  polymethineimine.  It  is  in  a  plane 
and  its  unit  cell  is  NCH.  The  five  independent  internal  coordinates  are  rN=c,  i~c-n, 
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Figure  4.1.  The  structure  of  polymethineimine  and  the  five  independent  internal 
coordinates:  rN=c,  tc-n,  tc-h,      and  (3. 


rC-H,  a,      Among  them,  r^=c,  fc-N,  and  a  determine  the  unit  cell  length  a  by 


Since  analytical  gradients  are  not  yet  available  in  PLH93  or  other  periodic  HF 
packages,  numerical  gradients  are  used  in  our  geometry  optimizations.    We  take 


lations.  To  be  sure  that  our  optimized  geometries  are  reliable,  we  have  checked  our 
results  by  doubling  these  step  sizes  and  reducing  them  by  half  to  assess  the  sensitivity 
of  the  numerical  results.  We  found  they  do  not  change  the  accuracy  of  our  reported 
data. 

Table  4.1  compares  our  optimized  HF  geometries  obtained  with  basis  sets  STO- 
3G  and  6-31G  for  polymethineimine  with  those  given  by  Teramae  et  al91  and  Karpfen.90 
From  the  table,  we  can  see  that  there  is  very  good  overall  agreement  between  our 


(4.22) 


ArN=c,  ArC-N,  ArC-//  to  be  0.01  A  and  Aa,  A/?  to  be  0.4  degree  in  our  calcu- 


6G 


Table  4.1.  Comparison  among  optimized  HF  geometries  with  basis  sets  STO-3G  and 
6-31G  for  polymethineimine  (Units:  A  and  Degree).  Our  alculation  was  done  by 
using  PLH93  with  21  unit  cells. 


TN=C 

rc-N 

rc-H 

a 

P 

STO-3G 

Ours 

1.281 

1.452 

1.097 

116.6 

125.0 

Teramae  et  al91 

1.280 

1.465 

1.107 

116.9 

125.3 

6-31G 

Ours 

1.266 

1.384 

1.088 

120.3 

121.8 

Karpfen90 

1.259 

1.388 

1.095 

121.0 

results  and  theirs.  But  we  also  discern  small  differences.  These  differences  may  come 
from  the  different  lattice  summation  cutoffs.  As  we  mentioned  above,  we  take  21  unit 
cells  in  our  calculations  plus  a  multipole  expansion  outside  this  region  to  enhance  the 
convergence  of  the  lattice  summations.  This  should  be  contrasted  with  only  9  unit 
cells  taken  in  Teramae  et  al's  calculation  and  five  in  Karpfen's. 

Table  4.2  lists  the  optimized  geometries  for  both  HF  and  MBPT(2)  methods 
with  three  basis  sets,  STO-3G,  6-31G,  and  6-31G**.  The  bond  length  rC-H  does 
not  change  much  with  basis  and  method  compared  to  the  other  internal  coordinates. 
With  increase  of  the  basis  set  size,  the  two  bond  lengths  r^=c  and  rc-N  always 
decrease  no  matter  whether  HF  or  MBPT(2)  is  used  in  the  calculation.  The  larger 
the  basis  the  stronger  the  two  bonds.  The  decrease  of  the  two  bond  lengths  also  means 
an  increase  of  the  repulsion  among  the  nitrogen  nuclei,  which  causes  an  increase  in 
a.  From  Figure  4.1,  one  can  see  that  increasing  a  means  decreasing  j3.  The  data 
in  Table  4.2  confirm  our  above  arguments.  We  can  also  see  from  Table  4.2  that  the 
single  bond  between  nitrogen  and  carbon  atoms  is  more  sensitive  to  the  size  of  the 
basis  than  that  of  the  double  bond.  The  result  is  that  the  difference  of  the  two  bond 
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Table  4.2.  The  optimized  geometries  for  polymethineimine  using  both  HF  and 
MBPT(2)  with  three  basis  sets,  STO-3G,  6-31G,  and  6-31G**  (Units:  A  and  De- 
gree). 


tn=c 

rc-N 

rc-H 

Q 

P 

STO-3G/HF 

1.2814 

1.4515 

1.0973 

116.63 

125.03 

STO-3G/MBPT(2) 

1.3367 

1.4836 

1.1190 

113.76 

126.33 

6-31G/HF 

1.2657 

1.3839 

1.0876 

120.30 

121.82 

6-3lG/MBPT(2) 

1.3039 

1.4098 

1.1088 

118.31 

122.96 

6-3lG**/HF 

1.2616 

1.3617 

1.0939 

119.61 

121.55 

6-31G**/MBPT(2) 

1.2859 

1.3717 

1.1023 

118.19 

122.16 

lengths  decreases  with  the  increase  in  the  basis  set  size.  Since  electron  correlation 
typically  diffuses  the  electron  cloud,  the  correlated  bond  lengths  tend  to  be  larger 
than  those  obtained  by  HF  with  the  same  basis  set.  Because  of  the  increase  of  the 
two  bond  lengths,  rN=c  and  rC-N,  &  decreases  and  then  /?  increases  with  inclusion  of 
electron  correlation.  The  double  bond  is  more  sensitive  to  electron  correlation.  This 
also  results  in  a  decrease  in  the  difference  between  the  two  bond  lengths,  r^=c  and 
rc-N- 

4.3.2    Band  gap 

Polymethineimine  was  found  to  be  a  semiconductor  in  the  early  seventies,  al- 
though its  band  gap  was  not  accurately  measured.54'55  Since  it  has  a  conjugate  struc- 
ture similar  to  polyacetylene,  polymethineimine  became  a  candidate  for  a  possible 
conducting  material  with  "doping".  Using  the  HF  method  with  a  small  STO-3G 
basis,  Teramae  et  al  predicted  that  polymethineimine  had  a  very  large  band  gap, 
suggesting  that  it  would  not  be  a  good  candidate  system.91  In  the  following,  we 
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Table  4.3.  Band  gaps  calculated  using  HF  and  MBPT(2)  with  three  basis  sets,  STO- 
3G,  6-31G,  and  6-31G**  (Unit:  eV). 

STO-3G/HF  -7.3126  2.7502  10.0628 

STO-3G/MBPT(2)  -6.1099  2.1524  8.2623 

6-31G/HF  -9.0705  -0.8250  8.2454 

6-31G/MBPT(2)  -7.5709  -2.5668  5.0041 

6-31G**/HF  -8.9331  -0.3900  8.5431 

6-31G**/MBPT(2)  -7.6166  -2.8340  4.7826 

will  show  that  by  augmenting  the  basis  set  and  including  electron  correlation,  the 
theoretical  band  gap  will  reduce  dramatically. 

The  second  column  in  Table  4.3  lists  the  HF  and  MBPT(2)  quasi-particle  orbital 
energies  for  highest  occupied  orbital,  eVtTnax,  the  absolute  values  of  which  are  the  first 
ionization  potentials  (IP)  predicted  by  the  corresponding  combination  of  method 
and  basis  set.  The  IP  increases  with  an  increase  in  basis  but  decreases  with  the 
inclusion  of  electron  correlation.  The  third  column  gives  the  the  HF  and  MBPT(2) 
quasi-particle  orbital  energies  for  the  lowest  unoccupied  orbital,  eC)mjn,  the  negative 
values  of  which  are  theoretical  electron  affinities  (EA).  The  larger  the  basis  set,  the 
stronger  the  electron  attachment  predicted  by  the  theories.  It  can  also  be  seen  that 
electron  correlation  plays  a  key  role  in  electron  affinities.  The  IP  and  EA  predicted 
by  MBPT(2)/6-3lG**  are  7.6166  and  2.8340  eV,  respectively. 

Table  4.3  lists  the  band  gaps,  Eg,  obtained  by  using  HF  and  MBPT(2)  with 
the  three  different  basis  sets.  The  optimized  geometry  with  the  same  method  and 
basis  is  used  in  each  calculation.  The  number  of  unit  cells  in  the  lattice  summation 
is  the  same  as  that  used  in  geometry  optimization,  namely  21.  HF  with  a  STO-3G 
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basis  predicts  that  the  band  gap  is  10.0628  eV,  which  is  much  too  large  to  match  the 
semiconducting  property  of  polymethineimine.  By  increasing  the  size  of  the  basis  set, 
HF  with  6-31G  reduces  the  value  by  1.8174  eV.  But  there  is  no  further  improvement 
when  the  larger  basis  set,  6-31G**,  is  used.  By  including  electron  correlation  through 
MBPT(2),  the  band  gap  is  further  reduced  by  1.8005  eV  for  STO-3G,  3.2413  eV 
for  6-31G,  and  3.7605  eV  for  6-3 1G**,  respectively.  The  larger  the  basis  set  the 
larger  the  contribution  of  electron  correlation  to  the  band  gap.  The  MBPT(2)  band 
gap  obtained  with  6-31G**  is  4.7826  eV,  which  is  about  0.7917  eV  larger  than  the 
corresponding  band  gap  of  polyacetylene.43  Considering  the  measured  band  gap  for 
polyacetylene  is  about  1.8  eV,  from  our  prior  MBPT(2)  study,  we  can  estimate  that 
the  real  band  gap  of  polymethineimine  is  around  2.6  eV.  Hence,  polymethineimine 
might  be  a  reasonable  candidate  for  "doping" . 

4.3.3    Vibrational  frequencies 

Three  fundamental  vibrational  frequencies  of  polymethineimine  were  observed 
in  the  early  seventies  by  Wohrle.91  Besides  the  semiconducting  property,  they  are 
the  only  measured  properties  for  polymethineimine.  Using  the  HF  method  with  ba- 
sis set  STO-3G,  Teramae  et  al91  calculated  the  fundamental  frequencies  for  polyme- 
thineimine. Teramae  et  al's  numerical  and  Wohrle's  measured  vibrational  frequencies 
are  listed  in  the  first  and  third  rows  in  Table  4.4.  The  theoretical  results  were  scaled 
by  a  factor  of  0.89.  With  the  scaling,  Teramae  et  al's  numerical  results  match  the 
observed  values  reasonably  well.  However,  reliable  theoretical  results  can  only  be  ob- 
tained with  a  basis  set  that  is  large  enough  to  get  the  results  sufficiently  converged, 
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Table  4.4.  Comparison  among  Teramae  et  al,  ours,  and  experimental  results.  The 
two  groups  of  the  numerical  vibrational  frequencies  were  obtained  using  HF  with 
STO-3G  and  scaled  by  a  factor  0.89  (unit:l/cm).  Our  calculation  employs  21  unit 
cells  plus  multipole  expansion  outside  the  region  compared  to  9  in  the  Teramae  et  al. 


C-N  str. 

C-H  def. 

C=N  str. 

C-H  str. 

Teramae  et  al91 

1152 

1351 

1669 

3077 

Oursa 

1176 

1334 

1685 

3109 

Experiment55 

1410 

1620 

3170 

and  with  the  inclusion  of  electron  correlation.  In  the  following,  we  will  calculate  the 
fundamental  vibrational  frequencies  using  both  HF  and  MBPT(2)  with  three  basis 
sets,  STO-3G,  6-31G,  and  6-31G**. 

The  structure  of  all-trans  polymethineimine  is  described  in  Figure  4.1  and  the 
internal  coordinates  optimized  with  the  combinations  of  two  methods  and  three  basis 
sets  are  listed  in  Table  4.2.  Since  the  stable  structure  is  a  plane,  all  trans  polyme- 
thineimine has  Cs  symmetry  besides  its  translational  invariance.  Cs  symmetry  has 
two  irreducible  representations,  A'  and  A".  Then  the  normal  modes  of  polyme- 
thineimine can  be  categorized  into  A'  and  A"  symmetry. 

The  unit  cell  for  all-trans  polymethineimine  is  NCH.  Then  there  are  nine  Carte- 
sian coordinates  in  each  cell.  We  denote  them  as  xN,  yN,  zN,  xc,  yc,  zc,  xH,  Vh, 
zh,  taking  the  z  axis  to  be  perpendicular  to  the  plane  of  polymethineimine.  Then 
xn,  Vn,  xc,  Uc,  %h,  Vh  are  A'  representations  of  Cs  symmetry  while  zN,  Zc,  zH  are 
A"  representations.  For  fundamental  vibrational  modes,  translational  invariance  is 
always  kept.  Then  only  one  among  rC-N  and  a  is  independent  since  the  unit  cell 
length  is  known  and  fixed.  Therefore,  there  are  only  four  internal  coordinates,  rN=c, 
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rc-N,  ai  and  0,  belonging  to  A'  and  there  are  four  A'  normal  modes  which  have 
non-zero  frequencies.  It  is  easy  to  see  that  for  the  A"  representation  there  is  only  one 
internal  coordinate  7,  the  angle  between  the  carbon-hydrogen  bond  and  the  plane 
of  the  polymer.  Then  there  is  one  normal  mode  with  non-zero  frequency  for  A".  In 
the  following,  we  only  consider  the  normal  modes  for  A'  representations  since  their 
frequencies  have  been  measured. 

Our  HF/STO-3G  vibrational  frequencies  calculated  with  the  geometry  given  by 
Teramae  et  al  are  given  in  the  second  row  of  Table  4.4.  To  compare  with  Teramae 
et  al's  results,  the  frequencies  listed  in  the  row  have  been  scaled  by  a  factor  0.89  as 
they  did  for  their  data.  The  two  calculations  agree  well.  In  fact,  since  the  numbers  of 
unit  cells  used  in  the  lattice  summations  of  the  two  calculations  were  not  the  same, 
this  causes  some  diffirences.  Considering  this,  we  can  see  that  the  error  induced  by 
the  numerical  gradient  and  hessian  are  sufficient  small  for  our  current  application. 

Table  4.5  provides  our  calculated  fundamental  vibrational  frequencies  for  poly- 
methineimine.  The  vibrational  frequencies  listed  in  each  row  are  calculated  by  the 
combination  of  method  and  basis  given  in  the  first  column  at  the  corresponding  op- 
timized structure  listed  in  Table  4.1.  The  numerical  vibrational  frequencies  given 
in  Table  4.5  have  not  been  scaled  by  any  factor.  From  the  table,  we  can  see  that 
both  basis  set  and  electron  correlation  have  a  large  influence  on  numerical  vibrational 
frequencies.  The  larger  the  basis  set  the  better  the  numerical  results  compared  to 
the  measured  data.  We  have  carefully  checked  the  influence  of  the  step  size  used  for 


72 


Table  4.5.  Fundamental  vibrational  frequencies  with  symmetry  A'  and  calculated 
with  different  methods  and  basis  sets  (unit:l/cm).  The  results  are  not  scaled. 


C-N  str. 

C-H  def. 

C=N  str. 

C-H  str. 

STO-3G/HF 

1305.94 

1517.19 

1866.66 

3576.59 

STO-3G/MBPT(2) 

1014.92 

1352.07 

1505.99 

3340.30 

6-31G/HF 

1257.19 

1527.82 

1696.89 

3176.84 

6-3lG/MBPT(2) 

1173.12 

1396.08 

1596.59 

2957.92 

6-3lG**/HF 

1140.57 

1379.49 

1607.04 

3135.55 

6-3lG**/MBPT(2) 

1009.77 

1265.55 

1582.76 

2996.73 

Experiment55 

1410 

1620 

3170 

numerical  gradient  and  hessian  calculations  on  our  vibrational  frequencies.  The  error 
induced  by  numerical  derivatives  should  be  around  10  cm-1. 

For  STO-3G,  the  HF  vibrational  frequencies  are  larger  than  the  experimental 
values.  MBPT(2)  improves  the  numerical  results  by  a  few  hundred  wave  numbers. 
When  the  6-31G  basis  is  used,  both  HF  and  MBPT(2)  provide  good  agreement  with 
experiment.  In  fact,  the  agreement  is  much  better  than  we  expected  or  that  the  meth- 
ods can  achieve  in  finite  systems.  If  we  look  further  at  the  data,  we  find  that  the  HF 
results  are  larger  while  MBPT(2)  values  are  smaller  than  the  observed  frequencies. 
When  large  basis  set  6-31 G**  is  used,  MBPT(2)  frequencies  are  even  smaller  than 
those  obtained  with  6-31G,  although  MBPT(2)/6-31G**  frequencies  match  the  mea- 
sured values  reasonably  well.  If  one  just  looks  at  the  differences  between  them,  their 
agreement  is  at  the  level  that  MBPT(2)  can  obtain  in  finite  systems.  The  difference  is 
that  in  finite  systems,  the  vibrational  frequencies  calculated  with  MBPT(2) /6-31G** 
are  almost  invariably  larger  than  those  observed12  while  here  the  situation  is  the 
opposite.  HF/6-31G**  vibrational  frequencies  match  the  measured  values  perfectly. 
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Needless  to  say,  the  agreement  between  HF/6-31G**  and  experiment  does  not 
mean  that  HF  gives  a  better  description  of  the  system,  but  instead  suggests  that 
the  calculation  benefits  from  error  cancellation  in  the  basis  and  level  of  correlation. 
To  find  out  the  reason  for  this  unusual  behavior,  we  have  calculated  the  vibrational 
frequencies  for  several  methylenimine  oligomers  using  MBPT(2)  with  both  basis  STO- 
3G  and  6-31G**.  Recognizing  that  oligomer  structures  are  different,  which  makes 
the  analysis  only  approximate,  we  find  the  vibrational  frequencies  associated  with  the 
bonds  of  interest  decrease  with  the  increase  in  the  oligomer  length,  tebdering  toward 
the  values  for  infinite  systems,100  and  suggesting  that  the  fundamental  vibrational 
frequencies  in  polymethineimine  should  be  smaller  than  those  in  an  oligomer.  How- 
ever, comparing  the  measured  vibrational  frequencies  for  methylenimine101  and  those 
listed  in  Table  4.5  for  polymethineimine,  we  find  that  the  vibrational  frequencies  for 
the  latter  are  larger  than  the  corresponding  frequencies  in  the  former,  especially  for 
the  C  —  H  stretch  and  C  —  H  deformation.  This  means  that  there  is  some  other 
mechanism  which  increases  the  frequencies  of  the  chain's  normal  modes.  In  fact,  the 
experimental  vibrational  frequencies  for  polymethineimine  are  measured  in  a  situ- 
ation where  numerous  chains  screw  together.  The  hydrogen  bonds  formed  among 
hydrogen  atoms  in  one  chain  and  nitrogen  atoms  in  other  chains  may  make  both  hy- 
drogen and  nitrogen  vibrations  more  difficult  and,  thus,  increase  the  corresponding 
vibrational  frequencies.  Further  experiments  are  needed  to  clarify  ths  issue. 
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4.4  Conclusions 

Using  both  HF  and  MBPT(2)  with  basis  sets  STO-3G,  6-31G,  and  6-31G**, 
we  have  optimized  the  geometries  for  polymethineimine.  Both  basis  set  and  electron 
correlation  have  a  strong  effect  on  the  equilibrium  structure  of  the  system.  Either  ba- 
sis set  improvements  or  including  electron  correlation  reduces  the  difference  between 
the  single  and  double  carbon-nitrogen  bond  lengths. 

For  the  band  gap,  the  larger  the  basis  set,  the  smaller  the  band  gap  in  both 
methods.  However,  the  band  gap  obtained  by  HF  with  6-31G**  is  still  8.54  eV. 
Electron  correlation  has  a  dramatic  effect  on  the  band  gap,  with  MBPT(2)  giving 
4.78  eV  for  the  band  gap. 

We  have  also  calculated  the  fundamental  vibrational  frequencies  for  polyme- 
thineimine using  HF  and  MBPT(2)  with  three  basis  sets.  The  theoretical  results 
match  the  measured  values  very  well  when  a  large  basis  is  used.  The  comparison 
among  theoretical  and  experimental  vibrational  frequencies  also  indicates  that  the 
interaction  among  different  chains  may  be  non-negligible. 


CHAPTER  5 

MANY-BODY  PERTURBATION  THEORY  FOR  QUASI-PARTICLE  ENERGIES 


5.1  Introduction 

Ionization  potentials  and  electron  affinities  of  atoms,  molecules  and  polymers 
are  important  subjects  in  quantum  chemistry.  Hartree-Fock  orbital  energies  provide 
the  first-order  approximation  to  them  according  to  Koopmans'  theorem.63  However, 
the  Hartree-Fock  approximation  fails  in  providing  accurate  numerical  data,102-103 
due  to  the  lack  of  orbital  relaxation  and  electron  correlation  contributions.  Various 
methods  have  been  developed  to  provide  the  two  corrections  and  these  methods  have 
been  applied  to  numerous  systems  in  the  last  three  decades.13-15,19'20'43,104-106  The 
methods  can  be  classified  as  indirect  and  direct  approaches.  The  indirect  approaches 
perform  separate  Hartree-Fock  and  correlated  calculations  on  the  neutral  and  ion 
systems  and  subtract  the  resulting  total  energies  to  obtain  a  theoretical  ionization 
energy  or  electron  affinity.104,107  These  types  of  approaches  have  the  disadvantage 
of  substracting  two  nearly  equal,  large  numbers  plus  extensive  computational  effort. 
However,  they  do  potentially  provide  the  "most  complete"  calculations  since  both 
components  are  treated  as  fully  as  possible.  Direct  methods  avoid  these  deficiencies 
by  formulations  that  incorporate  the  initial  and  final  states  simultaneously.  The 
latter  means  using  common  one-particle  orbitals,  which  facilitate  the  cancellation  of 
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common  terms  before  calculation.  However,  such  methods  frequently  compromise 
the  complete  treatment  of  relaxation  e.g. 

Direct  methods  can  be  further  classified  roughly  as  three  types,  propagator13-18 
or  Greens  function  methods  (GFM), 19-23  many-body  perturbation  theory 
(MBPT),37'41,43'105'108'109  and  recent  equation  of  motion  (EOM)  coupled-cluster  (CC)106-110"112 
approaches  or  earlier  CI  type  methods.113  Propagator  methods  are  built  upon  the  self- 
energy  £(£),  which  is  energy  dependent.  They  also  contain  two  coupled  terms  for 
the  evaluation  of  either  IP's  or  EA's.  However,  it  has  been  shown15  that  the  two 
terms  can  be  decoupled  and  the  energy  dependence  removed,  subject  to  a  coupled 
cluster  reference  state,15  that  leads  to  IP(or  EA)-EOM-CC.  This  result  suggests  that 
a  finite-order  MBPT  approach,  which  has  to  converge  to  IP-EOM-CC  in  the  limit, 
should  have  the  same  advantage.  Through  second-order  for  both  orbital  relaxation 
and  electron  correlation,  Pickup  and  Goscinski  showed64  that  second-order  MBPT 
provides  equivalent  results  for  ionization  energies  with  the  second-order  quasi-particle 
approximation,  provided  that  a  diagonal  approximation  is  made  in  the  self-energy, 
i.e. 


42)  =  W})-  (5-1) 

Born,  Kurtz  and  Ohrn114  extended  this  equivalence  to  third  order.  Paldus  and 
Cizek105  gave  a  general  discussion  about  the  similarity  between  GFM  and  MBPT 
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while  Hernandez  and  Langhoff115  pointed  out  that  the  equivalence, 

4n)  (5-2) 

does  not  hold  for  fourth  and  higher  orders. 

In  this  chapter,  we  will  derive  a  general  expression  for  any  order  MBPT  cor- 
rection with  any  given  zeroth-order  one-electron  orbitals  for  quasi-orbital  energies 
defined  as  ionization  potentials  for  occupied  orbitals  and  electron  affinities  for  unoc- 
cupied orbitals.  There  are  two  types  of  diagrams  for  IP's  and  EA's.  The  first  type  of 
diagrams  are  linked  diagrams  which  can  be  further  separated  into  two  sets.  One  set 
consists  of  the  diagrams  having  one  and  only  one  bubble  with  a  fixed  index.  They 
mainly  account  for  orbital  relaxation  effects.  The  other  diagrams  are  derived  from 
the  diagrams  for  the  total  energy.  In  each  of  those  diagrams,  there  is  one  and  only 
one  line  with  a  fixed  index.  The  second  type  of  diagrams  are  the  diagrams  which 
have  two  or  more  lines  with  fixed  indices.  Because  of  mutual  cancellations  among  the 
contributions  of  the  orbital  relaxation  and  electron  correlation,  only  remaining  the 
diagrams  are  those  which  fall  apart  by  switching  any  two  ending  vertices  of  the  lines 
that  have  fixed  indices.  These  diagrams  can  be  expressed  as  unlinked  diagrams,  each 
separated  part  of  which  is  a  diagram  of  the  first  type.  The  relationship  between  the 
first  type  of  diagrams  of  MBPT  corrections  to  quasi-particle  orbital  energies  and  the 
total  self  energy  is  also  established.  It  is  shown  explicitly  that  Eq.(5.2)  does  not  hold 
for  fourth  or  higher  orders.  It  is  even  not  valid  for  third-order  when  non-Hartree- 
Fock  (HF)  orbitals  are  used.  We  will  also  show  that  the  orbital  relaxation  has  a  finite 
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contribution  to  the  band  energy  for  infinite  periodic  systems  although  its  effect  on 
the  zeroth-order  orbital  energy  and  total  energy  per  unit  cell  is  infinitesimal.  This  is 
important  since  it  is  often  said  that  relaxation  effects  in  infinite  systems  scale  as  l/N 
with  the  number  of  unit  cells  N,  and  are  thought  to  be  negligible  in  calculations. 

The  arrangement  of  this  chapter  is  as  follows:  In  Section  5.2,  a  brief  summary 
of  the  MBPT  formulas  for  the  total  energy  will  be  given  and  then  the  expressions 
for  ionization  potentials,  electron  affinities,  and  for  quasi-particle  energies  are  de- 
rived. In  Section  5.3,  the  rules  for  drawing  and  interpreting  the  diagrams  are  ex- 
plained. MBPT(l),  MBPT(2)  and  MBPT(3)  are  used  as  examples  for  using  these 
rules.  Then  conclusions  are  drawn.  Appendix  B  presents  pertinent  scaling  properties 
of  the  MBPT  diagrams  in  extended  systems. 


5.2    Many  Body  Perturbation  Theory 


5.2.1    Total  energy 


The  Hamiltonian  for  a  system  can  be  written  as 


t        i>j  "y 


(5.3) 


where  ht  is  a  one-particle  interaction  including  the  kinetic  energy  and  Coulomb  in- 


teraction with  the  nuclei.  Let 


/  =  h  +  u 


(5.4) 
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where  u  is  an  arbitrarily  chosen  one-particle  operator,  {\p}}  are  one-electron  wave- 
functions,  and 

f\p)  =  EfM-  (5-5) 

Our  convention  on  indices  is  that  p,  q,r,  s,  ...  are  unspecified  orthogonal  spin  orbitals, 
while  i,  j,  k,  I  are  occupied  and  a,  b,  c,  d  unoccupied. 

Then  the  second-quantized  Hamiltonian  of  the  system  in  the  space  of  {\p)}  is10'12 

H  =  H0  +  Y,fpiP]Q  ~  E(plfil9>Pf9  +  7  E^H™)^^,  (5.6) 

p^q  P9  pqrs 

where 

^o  =  E4Vp  =  E/pp^-  (5-7) 

V  V 

Using  |0)  to  denote  the  zeroth-order  ground  state  of  the  system  as  the  reference 
state  (the  Fermi  vacuum),  the  Hamiltonian  can  be  written  in  the  normal  product 
form116"118 


H  =  <0|ff|0)  +  W,  (5.8) 


where 
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p^q^sr 


(5.10) 


r, 


Ei<P»IM)  -  (p|%),  ifp  =  9, 


(5.11) 


The  one-particle  interaction  Vp9  becomes  fpq  if  &  is  chosen  as  the  Fock  operator. 
From  perturbation  theory,  the  total  energy  E  of  the  system  can  be  expressed 


as 
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E=(0\H\0)+  £  E{m). 


m>2 


According  to  MBPT,10,12 


(5.12) 


£M  =  (OK^o^-Aj^lIO), 


(5.13) 


where 


i-|o><oi 


(5.14) 


the  number  of  W  operators  is  m  while  L  means  the  limitation  to  linked  diagrams 
only. 
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5.2.2    Ionization  potential 

For  the  positive  ion  system  in  which  an  electron  in  the  i0th  orbital  is  vertically 
ionized,  we  can  take 

*i0  =  «o|0>  (5.15) 

as  its  zeroth-order  wave  function.  Now  the  reference  state  is  |0io)  =  and  the 
Hamiltonian  is 

Hlo  =  {0lo\H\0lo)  +  Uio  +  Wio,  (5.16) 

where 

(0lo\H\0lo)  =  (0\H\0)  -  4°)  -  £>0||n0>  +  ftM*),  (5.17) 

i 

Uio  =  -  £(w||m)wio  [p*q]  ,  (5.18) 

PQ 


pq  ^  pqrs 


(5.19) 
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The  extra  one-particle  interaction  U{0  accounts  for  the  orbital  relaxation  effect.  The 
total  energy  Ei0  can  be  expressed  as 


£ioH0io|£|0to)+£4m).  (5.20) 

m>2 


Same  as  for  the  neutral  system,  E^  can  be  calculated  by 


4m)  =  (Oiol        +  Wl0)Ri?(Ui0  +  Wl0)...Ri°(Ul0  +  Wi0)]L  \0lo),  (5.21) 


where 


EM-*,- Ho  K  ' 


Since  |0j0)  =  i0 10) ,  the  hole  index  i'  in  the  diagrams  of  Eq.(5.21)  exclude  i0  while  i0 
is  included  in  the  particle  index  a'.  Now  we  change  the  Fermi  vacuum  from  |0jo)  to 
|0),  i.e.  the  hole  index  from  i'  to  i  which  denotes  the  occupied  orbitals  in  |0)  and 
the  particle  index  from  a'  to  a  which  denotes  the  unoccupied  orbitals  in  |0).  Then 
Eq.(5.21)  can  be  rewritten  as 

E-™^   =   E^  +  | all  mth  —  order  linked  diagrams  either  having 
at  least  one  one  -  particle  interaction  Uio  or  having 
at  least  one  line  index  fixed  as  z0j.  (5.23) 
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Figure  5.1.  Ui0  can  be  replaced  by  a  dotted  bubble  which  has  a  fixed  index  i0  and  is 
associated  with  an  extra  minus  sign,  e.g.  a  and  b  are  equivalent. 

We  use  dotted  lines  to  denote  the  lines  with  fixed  indices  i0.  It  is  easy  to  get 
the  following  conclusions  for  the  dotted  lines: 

i.  Each  dotted  hole  line  brings  an  extra  minus  sign.  Since  a  minus  sign  is  already 
associated  with  a  hole  line  in  an  antisymmetrized  diagram,  there  is  no  net  con- 
tribution to  the  phase  factor  of  the  diagram  from  a  dotted  hole  line.  A  dotted 
particle  line  has  no  contribution  to  the  phase  factor  either,  since  no  phase  has 
been  associated  with  a  particle  line; 

it.  From  Eq.(5.18)  and  according  to  the  rules  of  antisymmetrized  diagrams,59,119  we 
can  express  the  one-particle  interaction  Uio  by  a  bubble  with  a  dotted  hole  line, 
e.g.  the  line  index  is  fixed  as  i0  and  an  extra  minus  sign  is  introduced.  Figure 
5.1  provide  the  diagrammatical  description  of  this  equivalence; 
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Hi.  The  dotted  lines  associated  with  any  vertex  cannot  be  more  than  two,  otherwise 
the  value  of  the  diagram  vanishes,  since  its  cancellation  would  be  caused  by  the 
exchange  of  two  equivalent  orbitals; 

iv.  For  the  same  reason,  one  line  must  be  outgoing  while  the  other  is  in-going  if  two 

dotted  lines  would  be  associated  with  the  same  vertex; 

v.  Therefore,  any  two  dotted  lines  can  never  be  equivalent  in  a  non-zero  diagram; 

Exchanging  the  ending  vertices  for  two  dotted  lines  always  give  a  diagram,  which 
is  different  from  the  original  diagram.  The  latter  can  be  an  unlinked  diagram. 

vi.  A  dotted  line  should  always  be  counted  as  a  line  which  is  not  equivalent  to  any 

other  line  in  determining  the  numerical  factor  although  its  original  line  would 
be. 

vii.  For  a  diagram  in  Figure  5.2,  i0  should  not  be  included  in  the  index  of  the  hole 
line  since  the  diagram  in  Figure  5.3  diverges  and  the  diagram  is  excluded  in 
the  second  term  of  Eq.  (5.23).  This  happens  when  and  only  when  the  original 
diagram  has  a  structure  which  is  shown  in  Figure  5.4,  where  there  is  no  time 
overlap  between  the  upper  and  lower  parts.  In  fact,  i'  should  be  kept  if  a  diagram 
as  described  in  Figure  5.5  is  introduced  when  it  is  replaced  by  i. 

The  two  diagrams,  a  and  b,  in  Figure  5.6  are  identical  except  the  parts  which 
are  explicitly  drawn.  According  to  the  interpretation  rules  for  an  antisymmetrized 
diagram,  their  algebraic  expressions  are  identical  except  for  their  phase  factor.  Since 
diagram  a  has  one  more  closed  circle  than  that  of  diagram  b  and  considering  that 
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Figure  5.2.  The  diagram  is  composed  of  two  parts  which  have  no  time  overlap  and 
are  linked  only  by  a  dotted  line  with  fixed  index  i0  and  a  hole  line.  For  this  type  of 
diagram,  i0  is  excluded  from  the  index  of  the  hole  line. 


io  A      T  io 


r  -  T 


Figure  5.3.  Same  diagram  as  that  in  Figure  5.2  but  the  hole  line  is  replaced  by  a 
dotted  line  too.  This  diagram  diverges  since  the  denominator  between  the  two  parts 
vanishes. 
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Figure  5.4.  The  original  diagram  of  the  diagram  in  Figure  5.2. 
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Figure  5.5.  The  diagram  is  composed  of  two  parts  which  have  no  time  overlap  and 
are  linked  only  by  dotted  lines  with  fixed  index  i0  and  a  hole  line.  For  this  type  of 
diagram,  i0  is  excluded  from  the  index  of  the  hole  line. 
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there  are  no  contributions  to  the  phase  factor  from  dotted  lines,  the  two  diagrams 
have  different  signs  and,  then,  they  cancel  each  other.  With  the  same  argument,  we 
know  that  the  two  diagrams  in  Figure  5.7  cancel  each  other.  Now  let  us  consider  a 
general  case  where  the  two  diagrams  have  structures  as  shown  in  Figure  5.8.  The 
vertices  for  the  two  dotted  lines  are  not  necessarily  different  from  each  other.  Some 
of  them  may  superimpose  as  those  in  Figures  5.6  and  5.7.  The  upward  arrows  in 
the  two  diagrams  do  not  mean  they  must  be  particle  lines  .  Each  of  them  can  be 
either  a  hole  or  a  particle  line.  It  is  easy  to  see  that  the  algebraic  expressions  for 
the  diagrams  are  the  same  except  their  phase  factors.  Since  a  dotted  line  makes  no 
contribution  to  the  phase  factor,  the  difference  of  the  two  diagrams'  phase  factors 
are  totally  determined  by  the  difference  of  their  circle  numbers.  Since  switching  the 
ending  vertices  of  two  lines  always  either  bring  one  more  or  one  less  closed  circle, 
the  two  diagrams  have  different  phase  factors.  Then  the  two  diagrams  cancel  each 
other  even  if  one  of  them  is  unlinked.  From  the  above  argument,  we  can  conclude  for 
diagrams  which  have  two  or  more  dotted  lines  that 

viii.  If  two  diagrams  can  be  interchange  by  exchanging  the  ending  vertices  of  two 
dotted  lines,  they  cancel  each  other. 

Let  the  diagram  in  Figure  5.9  be  a  diagram  with  n  [n  >  2)  dotted  lines  from 
the  second  term  of  Eq.(5.23).  The  dotted  lines  are  not  necessarily  all  particle  lines 
as  drawn  in  the  diagram.  Some  of  the  starting  and  ending  vertices  of  the  dotted 
lines  may  superimpose.  Let  us  first  discuss  the  case  where  switching  the  ending 
vertices  of  any  two  dotted  lines  does  not  introduce  an  unlinked  diagram.  Then  the 
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Figure  5.6.  The  undrawn  parts  in  the  dotted  square  boxes  in  diagrams  a  and  b  are 
exactly  the  same.  The  two  diagrams  have  the  same  magnitude  but  different  signs. 
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Figure  5.7.  Diagrams  a  and  b  have  the  same  magnitude  but  different  signs.  The  three 
vertices  do  not  have  to  be  in  the  same  order  as  drawn  in  the  diagram. 
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Figure  5.8.  Diagrams  a  and  b  cancel  each  other.  The  two  dotted  lines  do  not  have 
to  be  particle  lines. 

diagrams  obtained  by  switching  the  ending  vertices  all  exist  in  the  second  term  of 
Eq.(5.23).  These  diagrams  have  the  same  magnitude  but  their  phase  differences  from 
the  original  diagram  are  determined  only  by  the  numbers  of  the  circles  formed  by  the 
lines  in  them.  By  all  possible  switchings  of  the  ending  vertices  for  the  n  dotted  lines, 
we  can  get  a  total  of  n!  different  diagrams  which  have  the  same  structure  except  the 
connections  of  the  dotted  lines.  Half  of  the  n!  diagrams  among  all  can  be  obtained  by 
switching  the  terminating  vertex  an  even  number  of  times,  while  the  other  half  can 
be  accomplished  by  an  odd  number  of  switchings.  Since  an  ending- vertex  switching 
for  any  two  dotted  lines  always  change  the  number  of  circles  by  one,  the  former 
n!/2  diagrams  have  the  same  phase  factor  as  that  of  the  original  diagram  while  the 
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Figure  5.9.  A  diagram  with  n  (n'  =  n)  dotted  lines  which  do  not  have  to  be  particle 
lines.  Some  of  the  starting  vertices  of  the  dotted  lines  may  superimpose  their  ending 
vertices. 

latter  n!/2  diagrams  have  an  opposite  sign.  Therefore,  all  these  n!  diagrams  mutually 
cancel. 

Now  we  consider  the  case  where  there  are  at  least  two  dotted  lines  in  the  diagram 
by  switching  the  ending  vertices  of  which  the  diagram  becomes  unlinked.  Then  if 
we  take  off  all  the  dotted  lines  in  the  diagram,  the  diagram  must  disintegrate  into 
two  or  more  pieces.  Otherwise,  all  interactions  in  the  diagram  are  already  connected 
by  non-dotted  lines  and  no  connection  pattern  of  the  dotted  lines  can  introduce  an 
unlinked  diagram.  Let  us  use  t  to  denote  the  number  of  the  separated  pieces  obtained 
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by  taking  off  all  dotted  lines  in  the  diagram.  In  each  piece,  there  is  at  least  one  ending 
vertex  and  the  same  number  of  starting  vertices  for  dotted  lines.  Some  of  the  starting 
and  ending  vertices  may  superimpose.  Let  n{  be  the  number  of  ending  or  starting 
vertices  in  the  zth  piece  and  the  pieces  are  numbered  in  such  a  way  that  nx  >  ...  >  nt, 
where  n  =  ni  +  ...  +  nt. 

Let  us  consider  the  case  where  nx  >  2  and  use  G  to  denote  the  enssemble  of 
all  linked  diagrams  obtained  by  connecting  the  t  pieces  with  dotted  lines.  These 
diagrams  all  exist  in  the  second  term  of  Eq.(5.23).  We  can  divide  the  diagrams  in  G 
into  groups.  In  each  group,  all  the  diagrams  have  the  following  common  characters: 
i.  The  number  of  dotted  lines  going  into  or  out  of  the  first  piece  is  the  same;  ii.  The 
connection  patterns  of  the  dotted  lines  outside  the  first  piece  are  all  the  same.  Then 
each  diagram  can  only  belong  to  one  group  and  there  is  no  overlap  among  different 
groups.  Let  us  concentrate  on  one  such  group.  Since  both  the  connection  pattern 
outside  the  first  piece  and  the  number  of  outgoing  and  in-going  dotted  lines  are  the 
same,  say  n\,  the  differences  among  these  diagrams  are  the  options  of  the  starting 
and  ending  vertices  for  the  dotted  lines,  which  connect  the  first  piece  and  outside, 
and  the  connection  patterns  of  the  remaining  dotted  lines  in  the  first  pieces.  It  is 
easy  to  see  that  the  number  of  the  diagrams  in  the  group  is  nx !ni!/(m  -  ni)!.  We 
can  further  divide  the  group  into  smaller  groups.  In  each  such  smaller  group,  the 
starting  vertices  of  the  n\  outgoing  dotted  lines  and  their  connection  pattern  outside 
of  the  first  piece  are  the  same  and  the  difference  among  the  diagrams  are  the  ending 
vertices  of  the  n\  in-going  dotted  lines  from  outside  and  the  m  -  n\  dotted  lines 
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Figure  5.10.  A  diagram  with  n  (n'  =  n,  n  =  2,3,4,...)  dotted  lines.  Without  the 
dotted  line,  the  diagrams  falls  apart  into  n  pieces.  These  pieces  are  linked  by  dotted 
lines  in  such  a  way  that  they  form  a  circle.  In  each  piece,  there  is  only  one  dotted 
line  coming  in  and  one  dotted  line  going  out.  No  dotted  line  starts  and  ends  in  the 
same  piece. 

connecting  the  vertices  inside  the  first  piece.  There  are  a  total  of  nx  dotted  lines  to 
choose  their  ending  vertices  and  there  are  n\\  number  of  choices.  Different  choices 
introduce  different  diagrams.  Then  there  are  a  total  of  diagrams  in  the  smaller 
group.  With  the  same  reason  used  above,  the  diagrams  in  the  smaller  group  mutually 
cancel.  Since  the  diagrams  in  each  smaller  group  mutually  cancel,  all  the  diagrams 
in  G  have  no  net  contribution  and  so  do  all  the  diagrams  with  ni  >  2. 

When  nx  =  1,  all  equal  one  and  t  =  n.  Then  the  diagram  has  a  structure 
shown  in  Figure  5.10.  One  cannot  get  a  linked  diagram  by  exchanging  ending  vertices 
of  the  two  dotted  lines  an  odd  number  of  times  but  could  with  an  even  number.  All 
such  diagrams  have  the  same  phase  and  do  not  mutually  cancel.  Then  all  diagrams 
that  remain  have  a  structure  as  in  Figure  5.10. 

From  the  above  discussion,  we  can  draw  the  conclusion  that 
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ix.  For  the  diagrams  which  have  two  or  more  dotted  lines,  the  only  remaining  struc- 
tures can  be  described  by  the  diagram  in  Figure  5.10  where  exchanging  of  any 
two  dotted  line's  ending  vertices  (or  starting  vertices)  introduces  an  unlinked 
diagram. 

The  diagrams  with  one  dotted  line  in  the  second  term  of  Eq.(5.23)  also  remain 
since  there  is  no  reason  for  them  to  cancel  each  other.  From  the  conclusions  obtained 
above,  Eq.(5.23)  can  be  rewritten  as 

_|_        mih  —  order  linked  diagrams  either  having  one 
bubble  with  a  dotted  line  and  a  fixed  index  io  or  having 
a  dotted  line  with  a  fixed  index  i0|  +  {all  mth  —  order  linked 
diagrams  having  a  structure  as  described  in  Figure  5.10  with  two 
or  more  dotted  lines  whose  indices  are  fixed  as  z0|,  (5.24) 

where  no  diagram  in  the  third  term  of  Eq.(5.24)  can  be  linked  when  the  ending 
vertices  of  any  two  dotted  lines  in  it  are  switched. 

The  ionization  potential  Iio  of  a  neutral  system  from  the  i0  orbital  is  the  energy 
difference  of  the  ion  system  and  the  neutral  system,  i.e. 

CO 

Ii0  =  El0-E=Yl  (4m)  -  £M)  •  (5.25) 

m=0 


T?  ("»)  _ 
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Ii0  can  be  expressed  as 


oo 


-(m) 
so 


(5.26) 


m=0 


Then  from  Eqs.(5.9),  (5.12),  (5.17),  (5.21)  and  (5.24),  we  get 


^(n0||n0)  +  (io\u\i0) 


(5.27) 


and  for  m  >  2 


Ii0      =    |aZi  ra£/i  —  order  linked  diagrams  either  having  one 

bubble  with  a  dotted  line  and  a  fixed  index  i0  or  having 
a  dotted  line  with  a  fixed  index  z0}  +  {a//  mth  -  order  linked 
diagrams  having  a  structure  as  described  in  Figure  5.10  with  two 
or  more  dotted  lines  whose  indices  are  fixed  as  iQ\.  (5.2* 


5.2.3    Electron  affinity 

For  the  negative  ion  system  in  which  an  electron  attaches  in  the  a0th  orbital, 
we  can  similarly  take 


*«o  =410) 


(5.29) 
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as  its  zeroth-order  wave  function.  After  a  similar  derivation,  one  can  get  the  following 
formulas  for  the  new  reference  state  \0ao)  =  fyao 

Hao  =  (Ojtf  |0ao>  +  Uao  +  Wao,  (5.30) 

where 

(Ojtf  |0ao)  =  <0|ff|0)  +  eg)  +  E(^ol^ao)  -  (a0\u\a0),  (5.31) 


Uao  =  E(aoPllao9)iVao  [pfg]  ,  (5.32) 

pq 


Wao  =  £  VPqNao  ftq]  +W  (PQ\\rs)Nao  [piqUf]  .  (5.33) 

pq  ^  pqrs 


Then  total  energy  Eao  can  be  expressed  as 


Eao  =  (0ao\Hao\0ao)  +  ^  (5.34) 

m>2 


where 


Ea?]  =  (0ao|  [(^o  +  Wao)Ra0°(Uao  +  Wao)...Rl°{Uao  +  <)]l  |0ao),  (5.35) 
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ft  =    )~^[.  (5.36) 


Just  as  for  the  ionization  potential,  E^  can  be  rewritten  as 


E^   =   E^  +  |  all  mth  —  order  linked  diagrams  either  having 
at  least  one  one  —  particle  interaction  Uao  or  having 
at  least  one  line  index  fixed  as  a0|  (5.37) 

by  changing  the  reference  from  a0|0)  to  |0).  Use  a  dashed  line  to  denote  the  lines 
with  fixed  index  a0.  It  is  easy  to  see  that  a  dashed  particle  line  has  been  associated 
with  a  minus  while  a  dashed  hole  line  is  not.  This  is  just  the  opposite  to  the  dotted 
lines  for  the  ionization  potential.  We  know  a  minus  sign  is  already  associated  with  a 
hole  line  while  there  is  no  phase  associated  to  a  particle  line  in  an  antisymmetrized 
diagram.59'119  Thus  a  dashed  line  in  a  diagram  always  brings  a  minus  sign  for  the  phase 
factor  of  the  diagram  no  matter  whether  it  is  a  dashed  hole  line  or  a  dashed  particle 
line.  From  Eq.(5.33)  and  according  to  the  interpretation  rules  for  antisymmetrized 
diagrams,120  the  one-particle  interaction  Uao  can  be  treated  as  a  bubble  with  a  dashed 
particle  line.  In  fact,  the  dashed  line  in  the  bubble  can  also  be  considered  to  be  a 
hole  line.  Whether  it  is  treated  as  a  dotted  particle  line  or  dotted  hole  line  does  not 
influence  our  following  discussion. 

For  a  diagram  with  n  (n  >  2)  dashed  lines,  say  the  diagram  in  Figure  5.11,  the 
dashed  lines  bring  a  phase  factor  (-l)n.  There  are  a  certain  number  of  diagrams 
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in  the  second  term  of  Eq.(5.37)  which  have  exactly  the  same  structure  with  that 
of  the  diagram  in  Figure  5.11  except  the  connection  patterns  of  the  dashed  lines. 
These  diagrams  have  the  same  magnitude  and  can  be  interchanged  by  switching  the 
ending  vertices  of  the  dashed  lines.  Since  the  phase  factor  (— l)n  associated  with  the 
n  dashed  lines  does  not  change  with  the  switchings,  phase  differences  among  these 
diagrams  are  totally  determined  by  differences  of  the  circle  numbers  formed  by  all  the 
lines  in  each  diagram.  We  can  have  exactly  the  same  discussion  about  these  diagrams 
as  we  did  for  the  diagrams  with  dotted  lines  above. 
Similarly,  Eq.(5.37)  can  be  rewritten  as 

E^    =  +  {a//  mth  —  order  linked  diagrams  either  having  one 

bubble  with  a  dashed  line  and  a  fixed  index  ao  or  having 
a  dashed  line  with  a  fixed  index  a0|  +  {a//  mth  —  order  linked 
diagrams  having  a  structure  as  described  in  Figure  5.12  with  two 
or  more  dashed  lines  whose  indices  fixed  as  a0|,  (5.38) 

where  the  structure  of  the  diagram  in  Figure  5.12  is  the  same  as  that  in  Figure  5.10 
except  all  the  dotted  lines  are  replaced  by  dashed  lines.  Using  the  same  reason  as 
argued  for  the  ionization  potential,  a0  should  be  excluded  from  the  summation  for 
the  particle  line  a  in  the  diagrams  given  in  Figure  5.13. 
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Figure  5.11.  A  diagram  with  n  (n'  =  n)  dashed  lines  which  do  not  have  to  be  particle 
lines.  Same  as  in  Figure  5.9,  some  of  the  starting  vertices  of  the  dashed  lines  may 
superimpose  their  ending  vertices. 


The  electron  affinity  Aao  of  an  electron  adding  into  the  a0  orbital  is  the  energy 
difference  of  the  neutral  system  and  the  negative  ion  system,  i.e. 


oo 

Aao  =  E  -  Eao  =  £  (£<■»>  -  Eg)) 

m 


(5.39) 
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Figure  5.12.  A  diagram  with  the  same  structure  as  that  of  the  diagram  in  Figure 
5.10  but  with  dotted  lines  being  replaced  by  dashed  lines. 
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Figure  5.13.  There  is  no  time  overlap  between  the  upper  and  lower  parts  in  both 
diagrams  a  and  b  which  are  linked  by  a  particle  line  and  a  dashed  line  or  several 
dashed  lines,  a'  ^  a0  should  be  excluded  from  the  index  of  the  particle  line. 
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Aao  can  be  expressed  as 


Aoc  =  £  4?-  (5.40) 


ao 

m=0 


Then  from  Eqs.(5.9),  (5.12),  (5.31),  (5.34)  and  (5.38),  we  get 

4?  =  -4^         <}  =  -  £(™o|K)  +  (a0\u\a0)  (5.41) 

i 

and  for  m  >  2 

AW   =   -{all  mth  —  order  linked  diagrams  either  having  one 

bubble  with  a  dashed  line  and  a  fixed  index  a0  or  having 
a  dashed  line  with  a  fixed  index  a0}  -  [all  mth  -  order  linked 
diagrams  having  a  structure  as  described  in  Figure  5.12  with 
two  or  more  dashed  lines  whose  indices  are  fixed  as  a0}.  (5.42) 

5.2.4    Quasi-particle  energy 

Beyond  the  zeroth-order  approximation,  one  can  still  define  the  quasi-particle 
energies  as 


ei  =  -Ii  =  E-  E, 


(5.43) 
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for  occupied  orbitals  and  as 


=  -Aa  =  Ea  -  E  (5.44) 


for  unoccupied  orbitals.  They  are  clearly  no  longer  one-electron  energies.  Now  using 
heavy  lines  without  arrows  to  replace  the  dotted  lines  in  the  diagrams  for  the  ioniza- 
tion potentials  and  dashed  lines  in  those  for  the  electron  affinities  and  associating  a 
minus  sign  with  each  heavy  line  without  regard  to  whether  it  is  originally  a  hole  line 
or  a  particle  line,  we  can  express  the  quasi-particle  energy  for  orbital  p  as 

oo 

Cp  =  4°)  +  e(1)  +  E  cm  (5  45) 

171  =  2 

where 

ep]  =  ~  (Pl«b)  (5.46) 

and  the  explicit  expressions  for  e{,m>  (m  >  2)  can  also  be  obtained  from  Eqs.  (5.24) 
and  (5.38).  According  to  the  definitions  for  dotted,  dashed,  and  heavy  lines,  there 
is  no  phase  factor  associated  when  a  dashed  line  is  replaced  by  a  heavy  line  while  a 
minus  sign  must  be  attached  when  a  dotted  line  is  substituted  by  a  heavy  line. 
When  p  is  an  unoccupied  orbital, 

epm)    =    {al1  mth  ~  order  linked  diagrams  either  having  one 

bubble  with  a  heavy  line  and  a  fixed  index  p  or  having 
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a  heavy  line  with  a  fixed  index  p\  +  I  all  mth  —  order  linked 
diagrams  having  a  structure  as  described  in  Figure  5.12  with  two 
or  more  dashed  lines  whose  indices  are  fixed  as  pj  (5-47) 

for  m  >  2  according  to  Eq.(5.38).  This  expression  is  also  valid  for  m  =  1  if  we 
compare  Eq.(5.47)  with  Eq.(5.46).  Then  we  can  express  ep  as 

ep  =  ep0)  +ap  +  7p,  (5.48) 

where 

op   =   {all  linked  diagrams  either  having  one  bubble 

with  a  heavy  line  and  a  fixed  index  p  or  having 

a  heavy  line  with  a  fixed  index  pj,  (5.49) 


7P   =   {all  linked  diagrams  having  a  structure  as 

described  in  Figure  5.12  with  two  or  more  heavy  lines 

whose  indices  are  fixed  as  pj.  (5.50) 


For  a  diagram  having  a  structure  as  described  in  Figure  5.12  and  with  n  heavy 
lines,  its  magnitude  is  equivalent  to  that  of  a  diagram  as  given  in  Figure  5.14,  where 
the  diagram  falls  into  n  unlinked  pieces  and  each  piece  is  a  linked  diagram  with  one 
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heavy  line.  But  one  cannot  separate  the  n  pieces  into  two  groups  which  have  no 
time  overlap.  Otherwise,  the  whole  diagram  diverges  since  there  would  be  a  zero 
denominator.  Let  7^n)  be  the  collection  of  all  the  diagrams  described  in  Figure  5.14 
with  n  pieces.  Then 

7^   =   |a//  unlinked  diagrams  which  have  n  pieces, 

each  of  which  is  a  diagram  of  ap|.  (5.51) 

If  we  reconnect  the  n  pieces  in  Figure  5.14  by  the  heavy  lines,  there  are  (n-1)!  possible 
ways  to  connect  them.  Each  connection  gives  a  linked  diagram  as  described  in  Figure 
5.12.  All  the  (n  -  1)!  so  constructed  diagrams  are  different  linked  diagrams.  This 
means  that  there  are  (n-1)  diagrams  in  7^n)  which  give  the  same  unlinked  diagram 
when  they  fall  apart  in  the  way  as  shown  in  Figure  5.14.  Considering  that  there 
is  a  phase  factor  (-1)"-1  between  a  diagram  in  Figure  5.12  and  its  corresponding 
unlinked  diagram  described  in  Figure  5.14,  jp  can  be  expressed  as 

lP=E(-1)n~1(n-l)\1^.  (5.52) 

n>2 

When  p  is  an  occupied  orbital,  we  can  still  express  ep  as  a  summation  of  the 
three  quantities,  ep°\  ap,  and  Due  to  the  cancellation  among  the  two  minus  signs, 
one  of  which  is  introduced  in  Eq.(5.43)  while  the  other  of  which  comes  from  the 
replacement  of  the  dotted  line  in  each  diagram  with  a  heavy  line,  Eq.(5.49)  is  valid 
for  occupied  orbitals.  For  a  diagram  in  Eq.(5.51),  a  phase  factor  (-1)"  is  introduced 
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Figure  5.14.  Diagram  a  is  a  diagram  described  in  Figure  5.12.  There  is  no  heavy  line 
inside  A  but  inside  B  when  diagram  a  has  more  than  two  heavy  lines.  In  diagram  b, 
the  ending  vertices  of  the  two  heavy  lines  linking  A  and  B  are  exchanged  and  then 
A  and  B  becomes  unlinked.  If  B  has  more  than  two  heavy  lines,  it  has  a  structure 
described  by  Figure  5.12. 


when  all  the  n  dotted  lines  are  replaced  by  heavy  lines.  Then  from  Eqs.  (5.28)  and 
(5.43),  we  obtain 

7P  =  sumn>2  (n  -  1 ) .  (5.53) 
Combining  Eqs.  (5.52)  and  (5.53),  we  can  write 


7P  = 


=  < 


E„>2 (n  -  1)  !tJ">,  if  p  occupied, 

En>2(-l)n_1(«  -  Wpn),       if  P  unoccupied. 


(5.54) 
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<7P  can  also  be  expressed  as  sums  of  different  order,  e.g. 


oo 

a 


£  (5.55) 


where 


ap"^    ~  m^  ~  order  linked  diagrams  either  having  one 

bubble  with  a  heavy  line  and  a  fixed  index  p  or  having 

a  heavy  line  with  a  fixed  index  pj.  (5.56) 

From  the  Appendix  B,  we  know  that  diagrams  with  a  bubble  have  finite  con- 
tributions to  quasi-particle  band  energies  in  extended  systems,  although  their  con- 
tribution to  the  total  energy  per  unit  cell  approaches  zero  when  the  number  of  cells, 
N,  goes  to  infinity.  As  these  are  "orbital  relaxation"  terms,  i.e.  they  correspond 
to  changes  in  the  orbitals  in  the  independent  particle  reference,  there  is  a  residual 
orbital  relaxation  contribution  to  energy  bands. 

5.3    Diagrams  for  Quasi-Particle  Energies 

5.3.1    General  rules 

According  to  the  formula  in  Eq.(5.56),  we  can  get  the  following  rules  to  draw 
the  diagrams  for  where  m  >  1  and  p  can  be  either  an  occupied  orbital  or  an 
unoccupied  orbital.  There  are  two  types  of  diagrams.  One  is  the  diagrams  having 
one  and  only  one  bubble  with  fixed  index  p  and  the  other  is  the  diagrams  having  no 
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bubble  but  with  one  line  index  fixed  as  p.  The  following  provide  the  rules  for  making 
all  the  diagrams  for  cr£m)  and  interpreting  them  for  er£m): 

i.  Draw  all  mth-order  antisymmetrized  diagrams  which  have  one  and  only  one  bubble 

which  is  drawn  with  a  heavy  line.  Fix  the  index  for  the  bubble  in  each  diagram 
as  p. 

ii.  Draw  all  mth-order  antisymmetrized  diagrams  for  the  mth-order  correction  to 

the  total  energy  £(m) .  For  each  such  diagram,  fix  one  line  index  as  p  each  time 
and  get  all  distinct  diagrams.  The  lines  with  fixed  indices  p  are  drawn  by  heavy 
lines.  These  diagrams  are  those  for  ej,m> .  Remember  here  the  line  with  fixed 
index  p  is  not  equivalent  to  any  other  line  in  the  diagram. 

Hi.  If  the  two  parts  of  a  diagram  have  no  time  overlap  and  are  connected  only  by  a 
particle  line  and  a  hole  line  such  as  those  in  Figures  5.4  and  5.13a,  the  indices 
for  the  lines  cannot  be  p  at  the  same  time.  In  other  words,  when  the  index  of 
one  line  is  fixed  as  p,  the  index  of  the  other  line  should  exclude  p. 

When  the  diagrams  for  ajm)  are  interpreted,  we  use  the  rules  of  antisymmetrized 
diagrams120  with  the  following  two  modifications  for  numerical  and  phase  factors: 

iv.  When  the  numerical  factor  of  a  given  diagram  is  determined,  the  heavy  line  is 
counted  as  a  line  different  from  all  others  in  the  diagram. 

v.  When  the  phase  factor  of  a  given  diagram  is  determined,  a  minus  sign  is  asso- 

ciated to  a  heavy  line,  e.g.  the  heavy  line  is  always  counted  as  a  hole  line  no 
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matter  whether  the  line  is  a  particle  or  hole  line  in  the  original  antisymmetrized 
diagram  and  no  matter  whether  p  is  an  occupied  or  unoccupied  orbital. 

The  mth-order  diagrams  in  7P  can  be  easily  obtained  since  the  diagrams  for  jp 
are  composed  of  the  diagrams  of  av  with  all  possible  combinations. 

5.3.2    Connection  with  self  energy  in  GFM 

It  is  easy  to  see  that  each  diagram  for  av  has  a  corresponding  diagram  in  the 
mth-order  total  self  energy  ET  in  GFM  and  vice  versa.  The  heavy  line  in  a  ap  diagram 
corresponds  to  the  outgoing  and  in-going  lines  in  the  corresponding  GFM  diagram. 
With  the  diagonal  approximation  at  E  =  epm\  the  two  diagrams  give  exactly  the 
same  expression,  e.g. 


The  total  self  energy  can  be  divided  into  two  parts,  the  irreducible  self  energy  £ 
in  the  Dyson  equation  and  reducible  self  energy  For  canonical  HF  orbitals,  T,R 
begins  to  arise  in  fourth-order.  Then  ST  equals  £  in  the  approximations  lower  than 


that  7P  has  no  contribution  for  corrections  lower  than  fourth-order.  Then  from  Eqs. 
(5.48)  and  (5.57),  we  have 


(5.57) 


fourth-order.  Since  e£x)  vanishes  for  canonical  HF  orbitals,  we  know  from  Eq.  (5.54) 


(5.58) 
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a — -x 

Figure  5.15.  The  diagram  for  e^K  The  cross  represents  the  one-particle  operator  Vpq 
which  vanishes  for  canonical  Hartree-Fock  orbitals. 

43)  =       (40)))  .  (5-59) 

which  are  exactly  the  equations  given  by  Pickup  and  Goscinski64  and  by  Born,  Kurtz 
and  Ohrn11'1  .  respectively.  Since  7;)  has  non-zero  contributions  for  orders  higher  than 
third,  ep  and  £  cannot  hold  in  a  simple  relationship  as  given  by  Eq.(5.2).  In  fact,  it 
does  not  hold  even  for  third  order  when  non-HF  orbitals  are  used.  However,  Eq.(5.2) 
can  be  replaced  by  Eq.(5.57)  which  holds  for  any  order. 

5.3.3    Diagrams  and  formulas  for  MBPT(l)  and  MBPT(2) 

According  to  general  rules  given  in  the  section  5.2A,  there  is  only  one  diagram 
cj,1'  which  is  given  in  Figure  5.15.  Its  algebraic  expression  is 

41*  =  Vpp  =  -  <Pl«b>-  (5-60) 

i 

For  canonical  orbitals,       vanishes  since  Vpp  does. 

For  tp2),  one  may  think  that  two  cr^  diagrams  can  form  an  unlinked  diagram 
in  7P.  However,  the  two  parts  in  such  a  diagram  cannot  have  time  overlap.  Then 
it  cannot  be  a  diagram  in  jp.  Therefore,  we  just  need  to  consider  diagrams  in  cr(2). 
There  are  two  types  of  diagrams  for  af\  one  set  with  a  bubble  and  the  other  without. 
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O  A-x 


'-— x       v  -O 


Figure  5.16.  Diagrams  having  one  bubble  for  ef\  This  represents  the  relaxation 
effect  and  does  not  approach  zero  in  extended  systems  if  non-HF  orbitals  are  used. 


Let  us  first  consider  the  diagrams  with  one  bubble.  We  can  only  get  two  such  di- 
agrams described  in  Figure  5.16.  According  to  Section  III.  A,  the  algebraic  expression 
for  them  is 


Via(ap\\ip)  (ip\\ap)Vt 


+ 


JO)  JO)  '  (0)  (0) 
tj         ta  t,-     —  to 


(5.61) 


These  two  diagrams  account  for  the  relaxation  effect.  They  vanish  when  canonical 
orbitals  are  used.  For  non-HF  orbitals,  their  contributions  do  not  approach  zero  in 
extended  systems  according  to  the  scaling  properties  shown  in  Appendix  B. 

Now  consider  diagrams  with  one  index  fixed  as  p.  There  two  antisymmetrized 
diagrams  are  described  in  Figure  5.17  for  the  total  energy.  From  them,  we  can  obtain 
the  four  diagrams  shown  in  Figure  5.18.  The  first  two  diagrams  have  one-particle 
interactions.  The  algebraic  expressions  give 


5p  =  E^-E^k  (5.62) 
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r— X 


X 


Figure  5.17.  Diagrams  for  second-order  MBPT  correction  to  the  total  energy,  e,g. 
for  E®. 


r  — X  a  X 


X         v— -X 


Figure  5.18.  Diagrams  obtained  from  those  for        with  one  line  index  fixed  for  ep2). 


where  the  rule  in  in  III.  A  has  been  used.  The  expression  for  the  other  two  diagrams 
in  Figure  5.18  is 

1  ^         |Mla6)|2  1  KpalMf 

'     2  £-  ejo)  +  £(o)  _  f(o)  _  £(o)     2a^.€(o)+e(o)_c(o)_e(o)-  (5-63) 

Then  for       we  have 


ep   —  Ap  +  Bp  +  Cp. 


(5.64) 


Ill 


When  the  canonical  Hartree-Fock  orbitals  are  used,  Ap  and  Bp  go  to  zero  and 


ep2^  can  be  written  as64 


(2)  _  1  \<pi\\ab)\*  1  \(pa\m> 

P    '  2  h  40)  +  ^  -  el0)  -  40)     2  53  40)  +  ^  -  e<°>  -  ef  ' 


After  summing  over  spin  indices,  Eq.(5.65)  becomes43 


F(2)        V  V  mPlUB)\*-Re[{PI\\AB){BA\PI)] 

p  M  ,    fo)      to)  (o) 


B 


i  ab  eP  +  er   —  eA  —  e 

,  ^  ^  2|(iM||JJ)|2  -  Re[(PA\\IJ){JI\PA)] 

rr*         f(0)  +  ^(0)   ^(0)  *(°>        '        [  ' 

u  a  tP  -t-  eA  —      —  ey 


where  P,  I,  J,  A  and  B  are  spatial  orbital  indices. 
5.3.4    Diagrams  of  MBPT(3) 

For  ef  \  we  need  to  consider  diagrams  in  erp3)  and  also  the  third  order  diagrams 
in  7P.  In  the  following,  we  just  give  the  concerned  diagrams.  It  is  easy  but  tedious 
to  get  their  algebraic  expressions. 

Let  us  first  consider  the  diagrams  for  of\  There  are  two  types  of  diagrams. 
The  first  type  have  one  bubble.  There  are  eighteen  such  diagrams  which  are  given  in 
Figure  5.19.  They  do  not  vanish  in  extended  systems.  Even  with  canonical  Hartree- 
Fock  orbitals,  the  first  six  diagrams  still  have  finite  contributions.  For  there 
are  fourteen  diagrams  drawn  in  Figure  5.20.  Among  them,  eleven  have  one-particle 
interactions.  For  the  first  twelve  diagrams,  we  can  get  four  diagrams  for  e£3)  from 
each  of  them.  For  the  last  two  diagrams,  each  introduces  three  diagrams  for  e(3). 
Then  there  are  in  total  forty  two  distinct  diagrams,  described  in  Figure  5.21,  with 
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Figure  5.19.  Diagrams  with  one  bubble  for  ep3\  There  are  eighteen  such  diagrams. 

one  particle-interaction  Vpq  and  one  heavy  line  for  ef\  From  the  other  three  diagrams 
in  Figure  5.20  with  only  two-particle  interactions,  we  can  get  12  diagrams  for  ep3K 
These  diagrams  are  given  in  Figure  5.22. 

Now  consider  the  third-order  diagrams  in  jp.  This  is  the  first  time  we  deal 
with  the  diagrams  in  -yp.  As  we  mentioned  in  Section  II.  D,  the  diagrams  in  7p  are 
unlinked  diagrams  composed  of  those  of  ov  in  all  possible  ways.  Figure  5.23  gives  all 
six  third-order  diagrams  in  yp.  These  diagrams  come  from  the  six  linked  diagrams, 
shown  in  Figure  5.24,  with  two  heavy  lines  and  a  structure  as  described  in  Figure 
5.12.  One  can  also  use  these  diagrams  as  examples  to  check  why  there  are  different 
signs  for  occupied  and  unoccupied  orbitals  in  Eq.(5.54). 

When  the  canonical  Hartree-Fock  orbitals  are  used,  all  the  diagrams  with  one- 
particle  interactions  Vpq  vanish.  It  is  easy  to  check  that  the  remaining  diagrams  have 
a  one-to-one  correspondence  with  those  for  third-order  GFM.20'106,120 
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Figure  5.20.  Diagrams  for  third  MBPT(2)  correction  to  the  total  energy,  e.g.  for 

5.4  Conclusions 

A  general  expression  for  any  given  order  MBPT  correction  to  the  quasi-particle 
energy  is  given.  The  diagrams  can  be  classified  into  two  types.  The  first  type  of 
diagrams  are  linked  diagrams  which  can  be  further  separated  into  two  sets.  One  set 
consists  of  the  diagrams  having  one  and  only  one  bubble  with  a  fixed  index.  They 
mainly  account  for  orbital  relaxation  effects.  The  diagrams  of  the  other  set  are  those 
derived  from  the  diagrams  for  the  total  energy.  In  those  diagrams,  there  is  one  and 
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only  one  line  with  a  fixed  index.  The  diagrams  of  the  second  type  can  be  expressed  as 
unlinked  diagrams,  each  unlinked  part  of  which  is  a  diagram  of  the  first  type.  Since 
each  unlinked  part  has  one  line  with  fixed  index,  these  unlinked  diagrams  are  size 
intensive. 

The  rules  to  interpret  the  diagrams  for  the  quasi-particle  energies  are  explained. 
They  are  similar  to  those  for  antisymmetrized  diagrams  but  with  two  modifications 
for  numerical  and  phase  factors.  The  diagrams  and  formulas  for  the  second-order 
correction  and  the  diagrams  for  the  third-order  correction  are  given  as  examples. 
One  can  easily  obtain  the  diagrams  for  higher  orders. 

The  orbital  relaxation  has  finite  contributions  to  quasi-particle  band  energies 
although  its  contribution  to  zeroth-order  orbital  energies  and  total  energy  per  unit 
cell  approaches  zero.  They  begin  to  have  a  contributions  in  second-order  if  non-HF 
orbitals  are  used. 

An  explicit  relationship  between  the  contributions  of  the  first  type  of  diagrams 
in  MBPT  corrections  to  quasi-particle  orbital  energies  and  the  total  self  energy  with 
the  diagonal  approximation  for  E  =  ej0)  in  propagator  or  Greens  function  methods 
is  also  established. 
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Figure  5.21.  Diagrams  without  bubbles  but  one-particle  interactions  for  e(3).  There 
are  forty-two  such  diagrams. 
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Figure  5.22.  Diagrams  without  bubbles  and  one-particle  interactions  for  e<3).  There 
are  twelve  such  diagrams. 
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Figure  5.23.  Diagrams  for  A^.  They  are  the  third-order  diagrams  in  7P.  There  are 
six  such  diagrams. 


Figure  5.24.  The  diagrams  of  with  two  heavy  lines  and  a  structure  described  in 
Figure  5.10  and  5.12.  These  diagrams  decompose  into  unliked  diagrams  described  in 
Figure  5.23. 


CHAPTER  6 

CONVERGENCE  OF  MANY-BODY  PERTURBATION  METHODS  WITH 
LATTICE  SUMMATIONS  IN  EXTENDED  SYSTEMS 

6.1  Introduction 

It  is  well-known  that  many-body  perturbation  methods  have  correct  scaling 
properties,59  size  extensivity121  and  size  intensivity,  its  analog  for  energy  differences, 
based  upon  the  linked  diagram  theorem.  However,  the  proof  was  given  for  two 
special  cases.  One  is  for  infinite  translationally  invariant  systems  such  as  nuclear 
matter122  or  the  homogeneous  electron  gas.4,123  The  other  is  the  supersystem  of  non- 
interacting  molecules,121  and  by  association,  the  interacting  ion.  Recently,  Deleuze  et 
al.  discussed  the  scaling  properties  for  oligomer  systems  of  increasing  size124  and  for 
inhomogeneous  one-dimensional  periodic  systems  like  polymers,58  and  showed  that 
there  is  a  possible  divergence  problem  for  high-order  MBGF.  They  showed  that  the 
integrand  for  intergration  over  k-space  could  diverge  logarithmically  with  lattice  sum- 
mations at  special  reciprocal  lattice  vectors  k  in  one-dimensional  periodic  systems, 
and  the  divergence  begans  to  occur  in  the  third-order  diagrams.58  A  further  analysis 
of  the  divergence  has  also  been  given  by  Nooijen  and  Bartlett.125 
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The  divergence  with  the  lattice  summations  can  disappear  after  the  integration 
over  k-space.  However,  there  is  as  yet  no  rigorous  mathematical  proof  of  the  con- 
vergence. Such  convergence  is  the  key  to  establishing  the  correct  scaling  properties 
of  many-body  perturbation  methods  in  inhomogeneous  extended  systems.  Besides 
the  theoretical  importance,  the  analysis  of  many-body  perturbation  methods'  conver- 
gence properties  with  lattice  summations  in  extended  systems  has  practical  impor- 
tance in  numerical  calculations,  particularly  with  regard  to  infinite-order  summations 
of  diagram  techniques  like  CC  theory.11  Since  it  is  impossible  to  do  an  infinite  sum- 
mation over  a  lattice  in  real  applications,  one  has  to  impose  a  cutoff  somehow.  The 
divergence  of  the  integrand  with  lattice  summations  at  some  special  k  values  implies 
that  meaningless  numerical  results  would  be  obtained  if  the  cutoff  is  implemented 
correctly. 

In  this  chapter,  we  present  a  rigorous,  general  proof  about  the  convergence  of 
any-order  MBPT,  CC  and  MBGF  with  lattice  summations  and,  therefore,  estab- 
lish the  scaling  properties  of  many-body  perturbation  methods  in  inhomogeneous 
extended  systems.  Our  proof  is  not  only  for  one-dimensional  systems  but  also  for 
two  and  three-dimensional  periodic  systems.  It  is  also  shown  that  in  actual  numer- 
ical applications,  one  should  take  a  universal  cutoff  for  all  lattice  summations  in  all 
the  diagrams  of  the  same  order  correlation  correction  and  ensure  that  the  value  of 
the  correction  has  converged  with  the  cutoff.  One  should  not  insist  upon  converged 
results  for  the  integrand  at  each  k-point,  as  that  could  lead  to  a  numerical  mistake. 
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The  chapter  is  organized  as  follows.  In  Section  6.2,  the  two-electron  integrals 
over  the  Bloch  orbitals  are  separated  into  two  parts.  One  of  them  has  a  short-range 
character  and  can  be  calculated  exactly.  The  other  has  a  long-range  character  and 
is  expressed  by  a  multipole  expansion.  In  Section  6.3,  we  will  show  the  convergence 
with  lattice  summations  for  any  given  diagram  of  MBPT  in  one-dimensional  periodic 
systems.  In  Section  6.4,  we  will  generalize  the  conclusions  obtained  in  Section  6.3  for 
one-dimensional  systems  to  three-dimensional  systems.  Section  6.5  will  discuss  the 
convergence  of  MBGF  and  CC  diagrams.  Conclusions  will  be  given  in  Section  6.6. 

6.2    Multipole  Expansion  of  the  Two-Electron  Integrals 

In  an  n-dimensional  (n  can  be  one,  two  or  three)  periodic  extended  system  with 
basic  translational  vectors  ai,  an,  the  one-electron  spatial  orbital  wave  functions 
are  Bloch  functions 

<Wr)    =   e8k»'rupkp(r),  (6.1) 

where  upkp(r)  is  translationally  invariant,  p  is  the  band  index  and  kp  is  the  reciprocal 
lattice.  A  reciprocal  lattice  k  can  be  expressed  as 

k   =   kibi  +  ...  +  knbn,  (6.2) 

where  (j=l,...,n)  are  determined  by  a,  •  b^  =  %  (i  =  1,  ...,n;j  =  1,  ...,n).  0pkp(r) 
can  be  further  expressed  as  the  summation  of  the  atomic  orbitals 
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=   ^EelkpRjE^kpXJaW,  (6.3) 


where 


xi(r)   =   Xa(r-Rj-Ra)  (6.4) 

are  atomic  orbitals  centered  at  RQ  in  the  jth  unit  cell  and  N  is  the  total  number  of 
the  unit  cells  approaching  +oo.  The  overlap  between  Bloch  orbitals  is  given  by32,43 

<<Wr)|<Wr)>   =   Skpkq^(C^yCpS^,  (6.5) 

where 

Ska0  =  £e<kNx°(r)|xj(r)>.  (6.6) 
The  two-electron  integrals  over  the  Bloch  functions  can  be  expressed  as43 

((pkp)(qkq)\(rkr){sks))    =   5kq,T{kr+k3_kp)G{pqrskpkrks)/N,  (6.7) 
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where  T  is  an  operator  which  moves  the  variable  x  back  to  the  first  Brillouin  zone  if 
it  is  out  of  the  region  and 

G(pqrskpkrks)   =      E    exp{  -  i[(kr  +  ks  -  k„)  •  R9  +  kr  •  R,.  +  ks  •  Rs]} 

x(A^(ri)A?kr+ka-kp)(r2  -  RJIftflA^r,  -  R,) 

xA^(r2-R,)) 

E    exp{i[kr-R1-(kr  +  ks-kp)-R2  +  ks-R3]} 
Ri,R2,R3 

a/3j0 

E    exp{i[kr  ■  Rx  -  (kr  -  kp)  R2  +  ks  •  R3]} 

Ri  ,R2,R3 

x  E  (c^p)*(cf  (kr+ks-kp))*c^cfs 

a/370 

x(a/3Rl|7R20R2+R3).  (6.8) 

It  is  well  known  that  the  individual  ((pkp)(qkq)\(rkr)(sks))  approaches  zero  for  ex- 
tended systems.  However,  G(pgrskpkrks)  has  a  finite  value  and  could  be  infinite  in 
special  cases.  The  possibility  of  G(pqrskpkrks)  having  an  infinite  value  is  exactly  the 
reason  why  we  discuss  the  divergence  problem  of  many-body  perturbation  methods 
here.  In  the  following,  when  we  talk  about  the  divergence  of  the  two-electron  inte- 
grals over  Bloch  orbitals  with  lattice  summations,  we  always  mean  the  divergence  of 
G(pqrskpkrks). 

Since  (a/?Rl|7R*0R2+R3)  exponentially  decreases  with  Rx  and  R3,  respectively, 
G(pqrskpkrks)  converges  very  fast  with  Ri  and  R3  and  there  is  no  convergence  prob- 
lem for  both  of  them.  It  is  safe  to  truncate  the  summations  over  them  in  numerical 
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calculations.  In  contrast  to  Ki  and  R3,  G(pqrskpkrks)  has  a  slow  convergence  with 
R2.  Let's  use  V  to  denote  the  assembly  of  all  the  R2  =  Ni&i  +  ...  +  A^a,,  which 
satisfy  the  condition  that  N,  <  N0  (i  =  l,...,n)  where  TVo  is  a  given  integer.  Then 
G(pqrskpkrks)  can  be  written  as 

G(pqrskpkrks)   =   G0(pqrskpkrks)  +  Gi{pqrskpkrks),  (6.9) 


where 


G0(pqrskpkrks) 


E    E  exp{i[kr-R,  -(kr-kp)-R2  +  k,-R3]} 

R2eVRi,R3 

X  E  (CakpY(Cf  (kr+k.-kp))»CrkrCJkJ 

a(3y0 

x(a{3Rl\jR26R2+K3).  (6.10) 


can  be  calculated  exactly  and 


Gi(pqrskpkrks)   =     E    E  exp{i(kr  ■  Rt  -  (kr  -  kp)  •  R2  +  k,  •  R3)] 

r2  evRi.Rs 

x  E  (cfp)*(cf (kr+k3-kp))*c^cfa 

a/370 

x(a/?Rl|7Ra0R2+R3)  (6.11) 


accounts  for  the  slow  convergence  of  G(pqrskpkrks)  with  lattice  summation  over  R2. 
The  symbols  e  and  6  mean  "belong  to"  and  "do  not  belong  to" ,  respectively. 
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The  two-electron  integrals  over  atomic  orbitals  have  the  following  form 

(a/?Rl|7R20R2+R3)   =    /  dv^xair,  -  Ra^fa  -  R/?  -  Rx),     1  . 

J  |r2  —  ri| 

xx7(r2  -  R7  -  R2)xe(r2  -  Re  -  R2  -  R3) 
=   ^  dr!rfr2Xa(ri  -  RQ)X/j(ri  -  R/j  -  Ri) 
xx7(r2  -  R2  -  R7)xe(r2  -  R2  -  R6  -  R3). 
|R2-[ri-(r2-R2)]|.  (6.12) 


When  N0  is  large  enough,  a  Taylor's  expansion  can  be  used  to  express  |R2  -  [ri  - 
(r.-R,)]!"1  as 


IR.-^Ar.-R.)]!  =  ^-r'-Vi  +  rY:VVi  +  ^ 


where 

^2  =   |R2|,  (6.14) 
r'   =  n-(r2-R2).  (6.15) 

Substituting  Eq.(6.13)  into  Eq.(6.12)  and  then  Eq.(6.12)  into  Eq.(6.11),  we  obtain 


G^pqrskrKK)   =   £%*r<%(*r+W.   £  ±-exp[i(k,  -  kr)  •  R2)] 

R2  IV 

ea;p[i(kp  -  kr)  •  R2)]  +  (pJpk^TCk.+k.-k^k 


i?2 
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_  jkpkr  yr(kr+ks-kp)k:3  _  jr(kr+ks-kp)ks  ikpkr 
upr     uqs  uqs  upr 


xexp[i(kp  -  kr)  ■  R2)]  +  (6.16) 


where  Sk^,  dklk2  and  Vk^  are  denned  as 


=   EiCt'YC^S^,  (6.17) 

a/3 


*k2   =  E(^kl)*^k2d^,  (6.18) 

a0 


^   =  EiCtrcf^lh  (6-19) 

Q/3 


with  the  dipole  momentum 


<ft»   =   E^kR(x^(r)|r|X^(r))  (6.20) 


R 


and  the  quadrupole  momentum 


=   £e*'R(X°(r)|rr|;#(r)).  (6.21) 

R 


^pg1  2  looks  very  similar  to  the  overlap  matrix  of  Bloch  functions.  In  fact,  it  is  the 
overlap  of  the  two  Bloch  functions  only  when  ki  =  k2.  When  kj  =  k2  =  k,  the 
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orthonormality  of  the  Bloch  function  gives 


SP7   =   <W  (6.22) 

It  is  easy  to  see  that  in  Eq.(6.15),  the  second,  the  third,  and  the  higher  terms 
converge  with  the  lattice  summation.  Then  the  convergence  of  the  two-electron  in- 
tegrals over  Bloch  orbitals  just  depends  on  the  first  term.  Therefore,  to  study  the 
convergence  of  the  two-electron  integrals,  we  can  approximate  Gi(pqrskpkrks)  as 

G^pqrskpKK)   =  <Spkr"k-^k-+k-k^kV(iVo,kp-kr),  (6.23) 


where 


/(JV0,k)   =    £  Lxp[ik-R)} 

R  6V 


=  ^2    y,  exp^ + "' +  kn^ 

jn>N0      j!>N0  blal  +  —3n&n\ 

9n  v       V*  cos  (kdi)-cos  (knjn) 
~        ^  '"'  ^   h-,ai  +    7  a  I  '  (6'24) 

jn>N0      j!>N0        Ulal  +  ■■■Jn^n\ 
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6.3    MBPT  in  Periodic  One  Dimensional  Systems 

6.3.1    Singularity  of  G(pqrskrk,ka)  and  convergence  of  MBPT(2) 

For  a  one-dimensional  periodic  system  with  basic  translational  vector  a,  the 
reciprocal  lattice  is  k  =  kb  where  b  is  determined  by  b  •  a  =  1  and  f(N0,  k)  becomes 


/(JV0,k)   =   -  ]T  -Lexp[ikj] 

a  \j\>N0  Ml 


2  1 

-  Yl  -cos{kj).  (6.25) 

a  j>N0  3 


Since 


£  cos(kj) 

j>N0 


=   l\Yl  exp[ikj]  +       exp[  ~  ikj\\ 


2 

3>N0  j>N0 
1 

< 


\sin(k/2)\ 


(6.26) 


and  the  series  Cj  =  1/j  approaches  zero  monotonically,  f(N0,k)  converges  accord- 
ing to  Dirichet's  Test126  except  at  k  =  0.  It  is  easy  to  see  that  /(iV0,  k)  goes  to 
infinity  when  k  =  0.  Therefore,  the  precondition  for  the  appearance  of  the  singu- 
larity of  G(pqrskpkrks)  is  that  there  is  no  momentum  exchange  during  the  interac- 
tion. Considering  the  orthonormality  of  the  Bloch  orbitals  described  by  Eq.(6.21), 
G(pqrskpkrks)  becomes  singular  when  and  only  when 

P  =  r,        kp  =  kr,        q  =  s,        kQ  =  ks.  (6.27) 
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a 
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c 


Figure  6.1.  Three  types  of  two-electron  interaction  vertices  whose  two-electron  inte- 
grals over  Bloch  orbitals  could  diverge. 

To  satisfy  Eq.(6.26),  the  diagram  for  the  two-electron  interaction  corresponding  to 
G(pqrskpkrks)  must  be  one  of  a,  b  and  c  in  Figure  6.1. 

Figure  6.2  gives  the  diagrams  of  second-order  many-body  perturbation  theory 
(MBPT(2))  for  the  total  energy  and  the  band  structure,  respectively.57  From  these 
diagrams,  one  can  see  that  there  are  no  such  two-electron  interaction  vertices  which 
are  the  same  as  those  described  in  Figure  6.1.  Thus,  all  G(pqrskpkrks)  converge  with 
the  lattice  summation  R2,  as  do  the  MBPT(2)  corrections  for  both  the  total  energy 
and  the  band  structure. 

6.3.2    Convergence  of  higher-order  MBPT 

Figure  6.3  gives  three  diagrams  of  the  third-order  MBPT  for  the  total  energy. 
It  is  easy  to  see  that  the  middle  vertices  of  the  three  diagrams  are  exactly  the  same 
as  those  of  a,  b  and  c  in  Figure  6.1,  respectively.  Therefore,  the  singularity  problem 
begins  to  occur  in  third-order  MBPT.58  More  than  one  two-electron  interaction  vertex 
could  have  the  forms  described  in  Figure  6.1  in  diagrams  for  higher-orders. 
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Figure  6.3.  Diagrams  of  MBPT(3)  for  the  total  energy.  The  corresponding  two- 
electron  integrals  of  the  middle  vertex  in  each  diagram  could  diverge. 
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Now  consider  a  general  case,  the  mth-order  (m  >  2)  MBPT  corrections  for 
both  the  total  energy  and  the  band  structure.  The  mth-order  MBPT  correction  for 
the  total  energy  is  given  by  linked  mth-order  antisymmetrized  diagrams,4'123  each  of 
which  is  a  sum  of  several  Gordstone  diagrams.59  The  mth-order  correction  to  the  band 
structure  is  given  by  the  linked  mth-order  antisymmetrized  (or  Goldstone)  diagrams 
with  one  line  index  fixed.57  In  the  following,  we  will  study  the  convergence  of  an 
mth-order  Gordstone  diagram  with  /  (/  >  0)  line  indices  fixed. 

The  algebraic  expression  of  a  mth-order  Goldstone  diagram  with  /  line  indices 
fixed  can  always  be  written  as 


where  W  is  the  volume  of  the  first  Brillouin  zone  (BZ),  £  sums  over  all  band  in- 
dices and  spin  variables,  D  denotes  the  denominator  consisting  of  zeroth-order  band 
energies  and  use  has  been  made  of  Eq.(6.7). 

D  has  no  singularity  for  the  total  energy  if  the  zeroth-order  band  gap  Ef^  is 
not  zero.  For  the  band  structure,  D  could  become  singular.43  However,  if  we  confine 
our  discussion  to  the  range 


/  = 


'£D(PMRiS1)...(PmQm\RmSm) 


(6.28) 


4? + 


(6.29) 
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where  e\J  and  e}J  denote  the  lowest  unoccupied  and  highest  occupied  zeroth-order 
orbital  energies,  respectively,  there  will  be  no  singularity  in  D  either.  Then  the 
singularity  of  the  integrand  only  arises  from  Gl  (i  =  1,  ..m). 

Let  ./V0  be  the  maximum  among  all  the  iVJ  which  defines  the  multipole  expansion 
range  for  Gj,  e.g. 

Gl   =  Gl0  +  G\ .  (6.30) 
Substituting  Eq.(6.29)  into  Eq.(6.27),  we  obtain 

m 

I  =  E/t,  (6.31) 
t=o 

where  It  is  the  sum  of  all  the  terms  in  which  the  number  of  the  G*i  is  t.  Then  I0  has 
a  finite  value  and  no  singularity  problem.  For  any  given  term  in  It  (t  >  0),  it  can 
always  be  written  as 

A  =   k^ThE  /B/km_m...  JBzdk1G[Gl..G\Gt+\..GZD 

=    {2n)m-M  E iykm.l+1..jyk1G\GlM\Gt0+l...G^D,  (6.32) 

where  k  =  kb  has  been  used. 
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First,  let  us  consider  the  case  where  /  =  0,  that  is,  no  line  has  a  fixed  index  in 
the  diagram.  Then  Eq.(6.31)  becomes 

A  =  p^E/^^-X^M"-^^1^  (6-33) 

Not  every  G\  has  a  singularity.  According  to  Eq.(6.26),  it  has  a  singularity  only 
when  the  corresponding  part  of  the  diagram  can  be  described  by  one  of  the  a,  b,  c 
interactions  in  Figure  6.1.  Let  t'  be  the  number  of  the  Gx  which  has  a  singularity 
(f  <  t).  Then  Eq.(6.32)  can  be  rewritten  as 

A   =    (2n)rn+i  E  £/W-  J^dhfiNo^-K,)... 

f(N0,kPt,  -kTt,)CD,  (6.34) 

where  C  is  the  remaining  part  of  G\G\...G\GtJ~lG1Q. 

If  the  two  reciprocal  lattices  kp  and  kr  of  any  one  f(N0,kp  -  kr)  in  the  above 
equation  are  forced  to  be  identical  by  the  conservation  of  the  reciprocal  lattices 
on  vertices,  the  diagram  diverges.125  However,  these  kind  of  digrams  always  appear 
in  pairs.  Their  divergent  parts  cancel  each  other  and  their  total  contributions  are 
converged  (see  Appendix  C).  In  the  following,  we  only  consider  the  cases  where  the 
two  reciprocal  lattices  in  any  /  function  in  Eq.(6.33)  are  not  forced  to  be  identical. 
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If  all  G\  having  singularities  happen  to  be  linked  together  such  as  a,  b  and  c  in 
Figure  6.4,  we  can  substitute  ki  by 


k'i   =   h-Ki,  (6.35) 


where 


Ki   =    £      +  (6-36) 


j=t+i 


into  Eq.(6.33)  and  then  get 


A  -  j2^^Ldkm+i- 


[  dkv+x 

J -7T 

rKi'  dk't,f(N0X)-  dk[f(N0,k[)CD.  (6.37) 

J-n-Kti  J-n-Ki 


Since  both  C  and  D  are  smooth  and  periodic  functions  of  ki,  A;m+1  with  a  peri- 
odicity 27r,  Eq.(6.36)  can  be  rewritten  as 

A  =  ^^E£r^i--£r^+i£r^/(^,i4')-" 

f  dk[f(NQ,k[)CD 

J  -7T 

=  jyKf(No,K)--- jyKHNoX)Z(k[^  (6.38) 


134 


k, 
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ii  ii  i 


a 


Figure  6.4.  Three  types  of  connections  when  all  the  singular  Gi  are  linked.  The 
directions  of  all  the  arrows  in  c  can  be  simultaneously  inverted. 


where 


Z{k[,  k'2, k'v)   =    (2;r)!L+i  E  //U-  ly^CD  (6.39) 


is  a  smooth  and  periodic  function  of  k[,      k[,  with  a  periodicity  of  2n. 
Consider 


B(k'2,...,k'tl)   =   ~jyk[f{N0X)Z{k\X,---X')-  (6.40) 


By  substituting  Eq.(6.24)  into  Eq.(6.39),  we  get 


B{k'2,...,k'v)  =  -T-M^-X), 

"  ]>N0  J 


(6.41) 
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where 


Cj(k2,...,k't,)   =   —  J\os(k[j)Z(k[,k'2,...,k'tl)dk[,  (6.42) 

are  the  Fourier  coefficients  of  Z(k[,  k'2, k'v)  with  respect  to  k[  in  the  range  of  [n,  n]. 
Since  Z(k[,  k'2, k'v)  is  a  continuous  function  of  k[,  we  find58,127 


Cj{k'2,...X>)   =   \vk.2  k[i[Z(k[,k'2,...,k't,)},  (6.43) 

J 


where  Vk^k<ti[Z(k[,k2,  ...,k't,)]  is  the  total  variation  of  Z(k[,k'2, k'v)  with  k[  over 
[-tt.tt]  at  given  k'2,  k'v.  Let  use  M  =  V[Z(k[,k'2, A;J,)]  to  denote  the  total 
variation  of  Z(k[,k'2, with  all  the  variables,  each  of  which  varies  in  the  range 
of  [-7r,  7r].  Since  Z(k[,  k'2, k't,)  is  a  continuous  function  with  all  the  variables  in  the 
range,  M  should  be  a  finite  number.  Now  we  have 


\cj(k'2,...,k't,)\<^-.  (6.44) 


Since 


-■  <  2^  (6.45) 

j>N0  J  j>N0  J 

and  the  left  side  of  Eq.(6.44)  converges,  we  can  conclude  that  the  summation  on  the 
left  side  of  Eq.(6.40)  not  only  converges  with  j  but  also  uniformly  converges  with 
j. 128  Thus,  B(k'2, k't.)  is  also  a  continuous  and  periodic  function  of  k'2,  ...k't,  with  a 
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periodicity  2n.  From  Eq.(6.44),  we  see  that  B(k'2, k't,)  converges  with  j  at  least  as 
fast  as  Remember  that  only  the  continuity  of  C  and  D  with  momentum  vectors 
has  been  used  above.  If  Z(k[,  k'2, k'v)  has  continuous  first  derivatives  with  respect 
to  k,  B(k'2, k't,)  would  converges  with  j  as  1/j2. 

Repeating  the  above  procedure  to  k'2,  k't,,  we  can  conclude  that  A  uniformly 
converges  with  all  lattice  summations  even  though  the  integrand  for  the  integration 
over  k  space  can  diverge  at  special  values. 

When  all  the  singular  Gx  do  not  directly  link  to  each  other  as  described  in  Figure 
6.4,  we  can  always  separate  them  into  several  groups  such  that  in  each  group,  all  G\ 
are  directly  linked.  Then  we  can  do  the  same  thing  for  each  group  as  we  have  done 
above  for  the  case  where  all  the  singular  GY  are  linked.  Then  in  this  case,  A  also 
uniformly  converges  with  all  the  lattice  summations. 

Since  each  A  uniformly  converges  with  the  lattice  summations,  In  uniformly 
converges  with  the  lattice  summations  and  so  does  any  mth-order  Goldstone  diagram. 
From  the  above  argument,  the  value  of  an  mth-order  Goldstone  diagram  with  no  fixed 
line  index  can  be  written  as 

I  =  NF,  (6.46) 

where  F  converges  with  lattice  summations.  The  diagrams  with  no  line  index  fixed 
are  those  for  the  total  energy.  Then  Eq.(6.45)  demonstrate  the  size-extensivity  of 
MBPT. 
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When  /  =  1,  e.g.  one  line  index  is  fixed  in  the  diagram, 


/  = 


1 


/TT  r 
dkm...  J 


TT 


dkiDGlG2...G 


in 


(6.47) 


(27T) 


m 


—H 


If  the  two-particle  integral  with  the  line  having  a  fixed  index  has  no  singularity, 
all  the  arguments  for  the  convergence  of  the  diagram  are  the  same  as  that  for  the 
case  where  /  =  0.  Otherwise,  the  situation  should  be  described  by  one  of  the  a,  b, 
c  interactions  in  Figure  6.5.  In  this  case,  we  further  separate  all  the  G\  functions 
into  two  groups  described  as  in  Figure  6.5.  Then  applying  the  above  techniques  to 
each  group,  we  can  show  that  A  converges  with  lattice  summations  in  both  groups. 
Therefore,  we  can  conclude  that  the  diagrams  with  one  line  having  a  fixed  index 
converge  with  lattice  summations  and  /  has  a  finite  value.  This  means  that  the 
diagrams  having  one  line  index  fixed  are  the  diagrams  which  are  size-intensive  for 
properties  such  as  for  ionization  potentials  and  electron  affinities. 

If  a  diagram  is  an  unlinked  diagram  as  described  in  Fig.  5.14,  the  algebraic 
expression  can  still  be  expressed  by  Eq.(6.47)  according  to  Eq.(B.12).  The  discussion 
for  linked  diagram  having  one  line  with  fixed  index  can  be  applied  for  each  unlinked 
piece  of  the  unlinked  diagram.  Then  the  diagrams  in  jp  are  size-intensive.  When 
/  =  2  or  larger,  we  have  the  case  that  the  two  momentum  vectors  can  have  the  same 
fixed  values  in  one  /  function.  Then  the  /  function  diverges  as  logN.  However, 
considering  the  factor  Nl~l  in  the  integral, 


/„  oc  lim  Nl~l  (logN)s   =  0, 


for     I  >  2, 


(6.48) 
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Figure  6.5.  Three  types  of  situations  where  the  line  having  a  fixed  index  connects 
the  two  singular  vertices. 
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where  s  (s  <  I)  is  the  number  of  the  /  functions  which  have  two  fixed  momentum 
vectors  with  the  same  value.  Eq.(6.47)  shows  that  for  any  diagram  having  two  or 
more  lines  with  fixed  indices  vanishes  when  N  goes  to  infinity. 

6.4    MBPT  In  Periodic  Three  Dimensional  Systems 

6.4.1    Singularity  of  G(pqrskrkrk<)  and  the  convergence  of  MBPT(2) 

For  three-dimensional  periodic  systems  (crystals), 


f(N0,k)  is  not  well  defined  since  the  second  term  in  Eq.(6.48)  oscillates  at  infinity. 
However,  the  contribution  of  this  oscillating  term  approaches  zero  as  1/R  after  inte- 
gration over  k.  Therefore,  we  do  not  need  to  worry  about  this  term.  With  this  in 
mind,  we  can  see  f(N0,k)  diverges  only  when  k  =  0.  Then  G{pqrskpkrks)  becomes 
singular  only  when 


/(iV0,k) 


ai  x  a2  •  a3A;  Jro 
4tt       1  r 


ax  x  a2  •  a3  k2 


[cos(kRo)  —  lim  cos(kR)]. 


(6.49) 


V  =  r,         kp  =  kr,         q  =  s,         kq  =  k. 


(6.50) 
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This  is  a  generalization  of  Eq.(6.26).  Diagrammatic  expressions  for  the  precondition 
of  the  singularity  of  G(pqrskpkrks)  for  three-dimensional  systems  are  the  same  as 
those,  described  in  Figure  6.1  for  one-dimensional  systems. 

Subject  to  the  same  reasons  as  given  for  one-dimensional  systems,  we  can  con- 
clude that  there  is  no  divergence  problem  in  MBPT(2)  for  three-dimensional  systems. 

6.4.2    Convergence  of  higher-order  MBPT 

There  is  no  much  difference  between  one-dimensional  systems  and  three-dimen- 
sional systems  when  the  convergence  of  MBPT  is  discussed.  We  can  employ  ex- 
actly the  same  procedure  used  for  one-dimensional  systems  in  Section  6.3  to  three- 
dimensional  systems.  However,  we  should  also  remember  that  the  lattice  vectors 
and  the  reciprocal  lattice  vectors  can  no  longer  be  described  by  scalars  but  by 
three-dimensional  vectors.  The  argument  from  Eq.(6.37)  through  Eq.(6.44)  for  one- 
dimensional  systems  should  be  modified  to  satisfy  three-dimensional  systems.  In  the 
following,  we  will  just  focus  on  this  difference. 

For  three-dimensional  systems,  we  can  follow  a  similar  discussion  to  that  for 
one-dimensional  systems  in  Section  6.3  to  finally  get 

A  =  ^Ej^/(*b,kt)...j^^  (6.51) 

where 

Z(k[,k'2,-X)   =    wZ*+i  E  fBzdkm+1...  JBzdkt,+1CD.  (6.52) 
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Defining 


BQ4,...X)  =  ^JBzdKf(N,K)z(K,K-X) 

/*7T  rn  pit 

=     (27r)3  J-x  ^  J-*  ^  I-    dklf(N'  klbl  +  /C2b2  +  *3bs) 

x    Z^bi +  A;2b2  +  A;3b3,k'2,...,k;,)  (6.53) 


and  substituting  Eq.(6.23)  into  Eq.(6.51),  we  obtain 


B(k'2,...x)  =  £  ■■      *  , 

R~iV  biai  +  J2a2  +  j3a3| 


where 


cjii2is(^2>  —  j^*1))  (6.54) 


rn  rn  rir 

(k'2,  ...,k't,)   =   —  J  dk3  j  dk2  J  dkxcos{kiji)cos{k2j2) 

^cos{k3j3)Z(klbl  +  k2b2  +  k3b3,  k'2, k't,),  (6.55) 


are  the  Fourier  coefficients  of  Z{klhl  +  k2b2  +  k3b3,  k'2, k'(,)  with  respect  to  ku 
k2,  k3.  It  is  easy  to  see  that 


\B(k'2,...,k[>: 


T   i  

Rjlv/  liiai  +j2a2  +  j3a3\ 


chhh  (^2 >  ••  •>  ) 


< 


(a1a2a3)1/3 


E  .1/3  1/3  .1/3CJU2J3(^2>  •••)kt') 
RjeV  Jl     32  33 


(6.56) 


With  the  same  procedure  as  used  in  Section  6.3  for  one  dimensional  systems,  one  can 
easily  show  that  B(k'2, k't,)  uniformly  converges126-128  with  ju  j2  and  j3  at  least 
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1/3  1 /3  1 /3 

as  fast  as  1/j/  ,  \/j2  and  ,  respectively,  and  so  does  B(k'2, k'f,).  In  other 
words,  B(k'2,  ...,k{/)  converges  with  the  radius  of  the  lattice  summation  range  R  at 
least  as  fast  as  l/i?1/3.  Because  of  this  uniform  convergence,  the  periodicity  and  con- 
tinuity of  Z(k\,  k2,...,  k'(,)  with  respect  to  k2,...,  k't,  is  transfered  to  B(  k'2, k't,) 
after  the  integration  over  k[  and  summation  over  the  lattice  summations.  By  re- 
peating the  above  procedure  to  k2,  k't,,  we  can  show  that  A  uniformly  converges 
with  all  lattice  summations.  All  the  arguments  after  Eq.(6.45)  for  one-dimensional 
systems  are  also  suitable  for  three-dimensional  systems. 

6.5    MBGF  and  CC  Methods 

The  convergence  of  MBFG's  diagrams  can  be  shown  in  the  same  way  as  that 
for  the  MBPT  diagrams  with  one  line  index  fixed,  which  are  the  diagrams  for  the 
ionization  potential  and  electron  affinity.  In  fact,  for  any  given  diagram  of  MBGF 
such  as  a  in  Figure  6.6,  there  is  always  an  MBPT  diagram  described  by  b  in  Figure 
6.6.  The  latter  is  exactly  the  same  as  the  former  except  the  line  with  fixed  index  has 
been  replaced  by  the  outgoing  and  ingoing  lines  of  the  former.  One  can  imagine  that 
the  outgoing  line  and  the  ingoing  line  in  MBGF's  diagrams  connect  at  infinity.  Then 
all  the  steps  used  for  b  can  be  followed  to  remove  the  divergences  of  a  if  there  are 
such  divergences.  Since  we  have  shown  in  Sections  6.3  and  6.4  that  any  mth-order 
(m  >  2)  MBPT  for  the  ionization  potential  and  electron  affinity  uniformly  converge 
with  lattice  summations,  no  matter  whether  one-dimensional  or  three-dimensional 
systems,  and  that  their  values  are  intensive,  we  can  conclude  that  any  mth-order 
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a  b 

Figure  6.6.  Connection  between  a  MBGF  diagram  (a)  and  its  corresponding  MBPT 
diagram  (b)  for  ionization  potential  or  electron  affinity,  which  has  a  line  with  a  fixed 
index  indicated  by  the  heavy  line.  The  remaining  parts  of  the  two  diagrams  in  the 
boxes  of  dashed  lines  are  exactly  the  same. 

MBGF  converges  with  all  the  lattice  summations  and  their  values  are  intensive  too, 
meaning  they  approach  finite  values  when  the  lattice  summations  go  to  infinity. 

In  CC  diagrams,  there  is  only  one  two-electron  interaction  which  can  possibly 
diverge  with  lattice  summations  when  its  reciprocal  lattice  vectors  satisfy  Eq.(6.26). 
In  some  diagrams,  the  reciprocal  lattice  vectors  of  the  Gi  are  forced  to  always  satisfy 
Eq.(6.26)  by  the  conservation  of  the  reciprocal  lattices  at  very  special  conditions.125 
When  this  happens,  these  divergent  diagrams  always  appear  in  pairs  and  their  di- 
vergent parts  cancel  each  other,125  eaxctly  the  same  as  those  of  the  MBPT  diagrams 
described  in  the  Appendix  C.  Since  there  is  only  one  two-electron  interaction  in  any 
CC  diagram,  the  CC  diagrams  can  be  considered  as  special  cases  of  MBPT  diagrams. 
Therefore,  CC  theory  uniformly  converges  with  the  lattice  summations  and  is  size 
extensive. 
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6.6  Conclusions 

By  using  the  multipole  expansion  for  the  two-electron  integrals  over  atomic 
orbitals,  we  have  shown  that  the  divergence  property  of  the  integrand  for  integration 
over  the  reciprocal  lattice  vector  k,  is  determined  by  the  first  term  in  the  expansion, 
the  effect  of  which  is  equivalent  to  the  potential  produced  by  a  charge  lattice  with 
phase,  no  matter  whether  the  system  is  one-  or  three-dimensional.  The  phase  is  a 
function  of  the  reciprocal  lattice  vectors.  It  becomes  zero  only  when  the  momentum 
exchange  is  zero.  With  zero  phase,  the  two-electron  integral  G  over  the  Bloch  orbitals 
diverges. 

Because  of  the  orthonormality  of  the  Bloch  orbitals,  the  singularity  of  two- 
electron  integrals  over  Bloch  orbitals  occurs  only  when  the  corresponding  ingoing 
line  and  outgoing  line  have  the  same  direction.  Since  there  is  no  such  two-electron 
interaction  in  the  diagrams  of  MBPT(2)  for  both  the  total  energy  and  the  band  struc- 
ture, MBPT(2)  has  no  singularity  problem  in  extended  systems  for  either  polymers 
or  crystals.  The  potential  singularity  problem  begins  to  occur  in  third-order  MBPT 
and  exists  in  higher-order  MBPT. 

The  diagrams  can  always  be  divided  into  three  types.  The  first  type  of  the 
diagrams  has  no  divergence  problem  at  all,  such  as  the  diagrams  of  MBPT(2).  The 
second  type  are  the  diagrams  which  really  diverge,  e.g.  the  value  of  each  goes  to 
infinity.  The  diagrams  which  belong  to  this  type  can  always  be  grouped  into  pairs 
and  the  divergent  parts  of  the  two  diagrams  in  each  pair  cancel  each  other.  The  third 
type  of  the  digrams  do  not  diverge  as  a  whole,  but  they  have  singular  integrands. 
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For  the  third  type  of  diagrams,  we  have  shown  that  each  diagram  uniformly 
converges  with  lattice  summations  after  the  integration  over  the  reciprocal  lattice 
vectors  k  because  of  the  cancellation  of  the  lange-range  contribution  with  different  k 
by  using  variable  replacements  and  the  asymptotic  behavior  of  the  Fourier  coefficients 
for  continuous  functions.  A  MBPT  diagram  of  the  third  type  converges  with  the 
radius  of  the  lattice  summation  range  R  at  least  as  fast  as  1  /Rl/n  where  n  is  the 
dimensionality  of  the  system.  In  fact,  it  could  be  much  faster  if  the  Bloch  orbitals 
and  zeroth-order  band  structure  have  continuous  first-  or  higher-order  derivatives 
with  respect  to  the  reciprocal  lattice  k. 

The  singularity  appearing  in  the  third  type  of  diagrams  also  appears  in  the 
remaining  parts  of  the  second  type  of  diagrams  after  the  cancellation  of  their  divergent 
parts.  Their  convergence  properties  with  lattice  summations  are  exactly  the  same  as 
those  of  the  third  type  of  diagrams. 

Our  proof  applies  equally  to  one-dimensional  and  three-dimensional  periodic 
systems.  Considering  that  two-dimensional  systems  are  special  cases  of  the  latter, 
the  conclusions  hold  for  two-dimensional  systems,  too. 

We  have  also  shown  that  the  diagrams  of  MBGF  and  CC  converge  with  lattice 
summations  after  integration  over  k  space. 

Our  proof  clearly  indicates  that  in  numerical  applications  of  many-body  pertur- 
bation methods  to  extended  systems,  one  should  not  simply  pursue  the  convergence 
of  the  integrand  with  lattice  summations,  but  the  convergence  of  the  correlation  cor- 
rection considered.  In  practice,  one  should  use  a  global  lattice  summation  cutoff  for 
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all  two-electron  integrals  over  Bloch  orbitals  in  all  the  diagrams  considered  and  be 
sure  that  the  numerical  value  of  the  correction  has  converged  with  the  cutoff. 


APPENDIX  A 

COMPARISON  AMONG  PERIODIC  HARTREE-FOCK  PROGRAMS 


Since  the  MBPT(2)  or  higher  order  corrections  depend  on  the  zeroth-order  wave- 
function  and  orbital  energies,  it  is  very  important  to  check  the  HF  results  before  the 
MBPT(2)  calculations.  In  this  Appendix,  we  compare  the  HF  results  obtained  by 
the  two  HF  packages,  which  we  can  access,  for  extended  systems  with  the  reported 
results  given  by  Suhai41.  Among  the  two  programs,  one  is  PLH93,  developed  by 
Andre's  group  and  only  for  one-dimensional  periodic  infinite  polymer  systems.32  The 
other  is  CRYSTAL92  which  is  developed  by  Dovesi  et  al  and  can  be  applied  for 
one-dimensional,  two-dimensional  and  three-dimensional  periodic  systems.31 

The  system  studied  by  Suhai  was  alternating  trans-polyacetylene.  The  atomic 
basis  set  was  STO-3G  and  the  geometry  used  is  Gl  in  Table  2.1. 

Table  A.l  lists  the  HF  total  energy  per  unit  cell  E^f  and  the  HF  band  gaps 
E?F  calculated  by  PLH92,  CRYSTAL92  and  reported  by  Suhai41.  In  the  calculations 
of  PLH93,  we  took  the  default  value  1(T7  for  the  cutoff  threshold  of  the  two-electron 
integrals  and  1(T8  for  the  SCF  convergence  threshold.  The  data  for  PLH93  listed  in 
Table  A.l  were  obtained  with  N  =  17,  where  N  is  the  number  of  unit  cells  used  in 
the  lattice  summations.  CRYSTAL92  has  a  more  complex  cutoff  criteria,  that  we  do 
not  describe  in  detail  here.  The  total  energies  of  one  unit  cell  obtained  by  the  three 
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programs  are  listed  in  the  second  column  of  Table  A.l.  PLH93  and  CRYSTAL92 
agree  to  the  first  six  significant  figures  while  they  agree  to  Suhai's  E^CF  to  one  less 
significant  figure.  The  third  column  of  Table  A.l  lists  the  HF  band  gaps  obtained 
by  the  three  programs.  The  difference  between  the  result  of  PLH93  and  that  of 
CRYSTAL92  is  about  0.01  ev  while  the  difference  between  PLH93's  and  Suhai's  is 
0.5  ev.  These  differences  between  Suhai's  results  and  those  of  PLH93  and  CRYSTAL 
cannot  be  diminished  by  adjusting  the  cutoff  of  the  lattice  summations  and  the 
thresholds  of  SCF  convergence. 

Table  A. 2  shows  the  convergence  of  PLH93's  E^f  with  the  number  of  unit  cells, 
N.  It  is  from  the  table  that  HF  results  can  be  considered  as  converged  values  when 
N=17  if  the  accuracy  is  not  higher  than  fifth  decimal. 
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Table  A.l.  Comparison  among  extended  Hartree-Fock  programs.  System:  alternat- 
ing trans-polyacetylene;  basis  set:  STO-3G;  structural  parameters:  Gl.  For  PLH93 
N=17. 


EHF  (au)  Eg  (ev) 

PLH93  -75.947633  8.72 

CRYSTAL92  -75.947643  8.71 
SUHAI  -75.947283  9.22 


Table  A.2.  Convergence  of  E%F  calculated  by  PLH93  with  number  of  unit  cells,  N. 
Geometry:  Gl;  Basis  set:  STO-3G. 


N=9  N=13  N=17  N=21 

E£F  (au)       -75.947622       -75.947631       -75.947633  -75.947634 


APPENDIX  B 

SCALING  PROPERTIES  FOR  ANTISYMMETRIZED  DIAGRAMS  IN 

EXTENDED  SYSTEMS 


In  extended  systems,  the  spatial  orbital  wavefunctions  are  Bloch  functions 

M'H-^E^'Aftr-R,),  (B.l) 

where  N  is  the  number  of  unit  cells  in  the  system.  The  two-particle  intergrals  can 
be  expressed  as43 

((n1k1)(n2k2)|(n3k3)(n4k4))  =  +k2,k3  +  k4)G(nln2n3n4k1k3ki),  (B.2) 

where  G^^nzn^k^)  converges  with  N  except  in  the  case  kx  =  k3.60  When 
ki  =  k3,  the  asymptote  of  G{nin2nzn4kikzk4)  with  N  is 

G,(n1n2n3«4kik3k4)  ~  logN.  (B.3) 

We  first  consider  an  mth-order  diagram  which  has  no  lines  with  fixed  indices 
and  two-particle  interactions  only.  In  such  a  diagram,  there  are  a  total  of  2m  lines. 
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Its  value  can  be  expressed  as 

I  =  J2C(P^Wris^iP^\\r2S2)-{pmqm\\rmsm),  (B.4) 

where  p{,  qu  n  and  s,  are  abbreviations  of  {njkjWj}  and  C  is  the  product  of  the 
denominator,  the  numerical  and  the  phase  factors.  Each  index  has  one  partner  such 
that  they  are  always  equal.  Thus,  the  number  of  independent  indices  is  2m  and  the 
algebraic  expression  of  the  diagram  looks  like 

E       C(PlP2||P3P4>....  (B.5) 

Pl,P2,—  ,P2m 

Since  each  term  in  the  sum  of  the  diagram  is  a  product  of  m  two-electron 
integrals,  there  are  a  total  of  m  delta  functions  to  confine  the  freedom  of  the  2m 
k  vectors.  Each  k  appears  in  two  and  only  two  of  the  m  delta  functions.  The  delta 
function  5(ki  +  k2,  k3  +  k4)  demands  that 

k!+k2  =  k3+k4.  (B.6) 

The  m  delta  functions  provide  m  such  equations  for  the  2m  k  vectors.  If  the  diagram 
has  s  unlinked  parts,  all  the  equations  like  Eq.(B.6)  can  be  grouped  into  s  groups. 
There  are  no  common  k  vectors  among  different  groups.  Each  k  vector  appears  once 
on  the  left  side  of  one  equation  and  appears  also  once  on  the  right  side  of  another 
equation  within  the  same  group.  Then  if  we  sum  all  the  equations  in  a  group,  we 
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find  the  two  sides  always  equal.  This  indicates  that  there  is  one  equation  which  is 
not  independent  in  each  group.  Then  there  are  m  —  s  independent  delta  functions  to 
confine  the  freedom  of  the  k  vectors.  These  delta  functions  reduce  the  number  of  k 
vectors  as  summation  indices  to  bem  +  s  from  the  original  2m.  The  summation  over 
each  k  vector  £k  can  be  replaced  by  the  integration  ^  J  dk,  where  W  is  the  volume 
of  the  first  Brillouin  zone.  Then  Eq.(B.5)  can  be  written  as 

/  =  NSF,  (B.7) 

where 

F  =  777^7     £     /  CG1G2...Gmdk1dk2...dkm+s  (B.8) 

m...n2m-l 

and  Eq.(B.2)  has  been  used.  Although  each  G  function  can  approach  infinity  at 
special  points  in  k  space  as  logN  when  N  -4  oo,  it  can  be  shown  that  F  converges 
with  N.60 

For  one-particle  integrals,  we  can  get 


((nk)H(n'k'))    =  5(k,k')((nk)|u|(nk)> 

=   S(k,k')F(nn'k),  (B.9) 
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where 

F(nn'k)  =  5:<AS(r)|u|AS(r  -  R,)>.  (B.10) 

From  Eq.(B.9),  we  know  that  adding  a  one-particle  interaction  in  a  diagram  provides 
no  new  independent  k  vector  and  then  it  does  not  change  the  scale  property  of 
the  diagram.  Then  the  scaling  property  of  a  diagram,  whether  it  has  one-particle 
interactions  or  not,  are  totally  determined  by  the  number  of  the  unlinked  parts  of 
the  diagram,  e.g.  Eq.(B.7)  is  also  valid  for  diagrams  having  one-particle  interactions. 

From  Eq.(B.7),  we  can  conclude  that  the  value  of  a  diagram  without  fixed 
indices  is  proportional  to  Ns.  For  all  linked  diagrams,  s  =  1  and  then  their  values 
are  proportional  to  the  size  of  the  systems.  This  type  of  diagrams  are  those  for  total 
energy  and  the  proportionality  of  their  values  to  system  size  is  just  the  manifestation 
of  the  energy's  extensivity.  For  unlinked  diagrams,  s  >  2  and  the  proportionality 
does  not  hold  any  more.  They  do  not  satisfy  the  extensivity  requirement  as  diagrams 
for  the  energy. 

Now  we  consider  a  diagram  which  is  linked  but  has  /  (/  >  1)  lines  with  fixed 
indices.  Let  l0  (l0  <  I)  be  the  number  of  the  pieces  which  the  diagram  would  fall 
apart  if  all  the  I  lines  with  fixed  indices  were  taken  away  from  the  diagram.  Since  the 
lQ  pieces  do  not  mutually  exchange  momentum.  They  behave  like  unlinked  diagrams 
with  l0  unlinked  pieces.  Then  /0  delta  functions  will  become  redundant.  Then  the 
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number  of  the  independent  k  vectors  is  m  +  Iq  —  I  and  we  can  get 


/  =  Nlo'lF. 


(B.H) 


For  the  diagrams  which  have  either  a  bubble  or  one  line  with  fixed  indices,  l0  =  I  = 
1.  For  the  diagrams  which  have  more  than  one  line  with  fixed  indices,  but  their 
structures  can  be  described  by  Figures  5.10  and  5.12,  /0  =  /.  For  these  types  of 
diagrams,  we  have 


This  shows  the  intensivity  of  the  diagrams.  These  diagrams  are  for  IP's  and  EA's 
which  are  intensive  properties. 

When  l0  <  I,  I  approaches  zero  when  N  ->  oo.  This  means  that  if  a  diagram 
is  still  linked  after  exchanging  the  ending  vertices  of  its  two  lines  whose  indices  are 
fixed,  the  value  of  the  diagram  approaches  zero  when  N  goes  to  infinity. 


I  —  F. 


(B.12) 


APPENDIX  C 
CANCELLATION  AMONG  DIVERGENT  DIAGRAMS 


There  are  situations  that  kj  and  k2  are  forced  to  be  identical  such  as  the  case 
showed  in  Figure  C.l.  In  these  cases,  the  line  indices  at  each  side  of  the  two  vertices 
denoted  as  V  in  both  diagrams  a  and  b  in  Figure  C.l  must  be  identical  according  to 
Eqs.(6.21)  and  (6.22).  We  also  know  from  Eqs.(6.21)  and  (6.22)  that  functions  Gx 
for  the  two  vertices  are  identical  when  N0  is  given. 

Among  mth-order  (m  >  2)  diagrams,  if  there  is  a  diagram  which  can  be  de- 
scribed by  a  in  Figure  C.l,  there  must  be  another  diagram  which  is  described  by  b 
with  exactly  the  same  remaining  part  of  a. 

If  the  remaining  parts  of  diagrams  a  and  b  are  identical,  the  algebraic  expressions 
for  the  two  diagrams  are  the  same  except  for  opposite  signs  and  the  G  functions  for 
the  two  V  vertices.  Since  Gx  functions  for  the  two  vertices  are  exactly  the  same  with 
any  given  N0,  the  divergent  parts,  induced  by  V,  of  the  two  diagrams  cancel  each 
other  because  of  their  opposite  signs. 

The  G0  functions  for  the  two  vertices  are  different,  Therefore,  the  total  contri- 
bution of  the  two  diagrams  are  usually  not  zero. 
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Figure  C.l.  A  pair  of  divergent  diagrams  (a)  and  (b),  the  divergent  parts  of  the  two 
diagrams  cancel  each  other.  The  remain  parts  in  the  two  boxes  defined  by  the  heavy 
dashed  lines  are  exactly  the  same. 

For  CC  diagrams,  the  situation  could  be  a  little  more  complex  because  the  part 
which  forces  lq  and  k2  to  be  identical  can  have  several  different  forms.  However,  the 
conclusion  for  divergent  diagrams  is  the  same  as  that  for  MBPT  divergent  diagrams. 
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